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13-  ffigl&cond  year,  we  pursued  the  development  of  a  reaction  zone  model  which 

can  be  used  in  the  numerical  simulation  of  turbulent  combustion.  We  focused  our  attention  on  the 
high  Reynolds  number,  high  Damkohler  number  case' in  which  the  reaction  zone  thickness  is  much 
smaller  than  the  thickness  of  the  small  scales  of  turbulence.  Furthermore,  we  treated  the  case  in 
which  the  entire  flame  thickness,  including  both  the  convection-diffusion  zone  and  the  reaction  zone, 
is  smaller  or  equal  to  the  small  scale  of  turbulence.  Accordingly,  we  developed  on  unsteady 
strained,  multistep  chemistry  and  diffusion  thin  flame  model  which  is  based  on  the  fundamental 
equations  governing  the  convection-diffusion-reaction  balance  within  the  flame  structure,  while 
assuming  that  the  strain  rate  remains  constant  within  the  flame.  Moreover,  we  compared  the  results 
of  a  computation  of  an  exothermic,  finite  rate  kinetics  model  of  a  reacting  shear  layer  with  those 
obtained  using  an  infinite  rate  kinetics  case.  This  comparison  established  the  fact  that  the  impact  of 
heat  release  on  the  flow,  mixing  and  dynamics  can  be  obtained  using  a  thin  flame  model.  In  the 
coming  year,  we  will  focus  on  linking  the  thin  flame,  finite  rate  chemistry  model  to  the  flow 
simulation  code  and  develop  the  proper  set  of  boundary  conditions  of  the  flame  model  in  order  to 
expand  its  regime  of  applicability. 
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SUMMARY 

During  the  course  of  the  second  year,  we  pursued  the  development  of  a  reaction  zone 
model  which  can  be  used  in  the  numerical  simulation  of  turbulent  combustion.  We  focused  our 
attention  on  the  high  Reynolds  number,  high  Damkohler  number  case  in  which  the  reaction  zone 
thickness  is  much  smaller  than  the  thickness  of  the  small  scales  of  turbulence.  Furthermore,  we 
treated  the  case  in  which  the  entire  flame  thickness,  including  both  the  convection-diffusion  zone 
and  the  reaction  zone,  is  smaller  or  equal  to  the  small  scale  of  turbulence.  Accordingly,  we 
developed  an  unsteady  strained,  multistep  chemistry  and  diffusion  thin  flame  model  which  is 
based  on  the  fundamental  equations  governing  the  convection-diffusion-reaction  balance  within 
the  flame  structure,  while  assuming  that  the  strain  rate  remains  constant  within  the  flame. 
Moreover,  we  compared  the  results  of  a  computation  of  an  exothermic,  finite  rate  kinetics  model 
of  a  reacting  shear  layer  with  those  obtained  using  an  infinite  rate  kinetics  case.  This  comparison 
established  the  fact  that  the  impact  of  heat  release  on  the  flow,  mixing  and  dynamics  can  be 
obtained  using  a  thin  flame  model.  In  the  coming  year,  we  will  focus  on  linking  the  thin  flame, 
finite  rate  chemistry  model  to  the  flow  simulation  code  and  develop  the  proper  set  of  boundary 
conditions  of  the  flame  model  in  order  to  expand  its  regime  of  applicability. 
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I.  OBJECTIVES 

The  objectives  of  this  work  are: 

(1)  To  develop  a  computational  framework  for  the  simulation  of  turbulent  combustion 
phenomena  by  deriving  equations  which  govern  the  flow  and  combustion  in  different, 
physically-distinct  regions  in  the  domain.  This  formulation  is  independent  of  the  numerical 
method,  or  methods,  used  to  integrate  these  equations  and  hence  will  have  applications  beyond 
vortex  simulations; 

(2)  To  develop  fundamentally  based  flame  structure  models  using  the  equations  obtained 
in  (1)  where  flow-combustion  interactions  are  represented  by  the  effect  of  the  time-dependent 
strain  rate  exerted  by  the  outer  flow  on  the  flame  structure.  One  important  aspect  of  these 
models  will  be  the  incorporation  of  multistep  chemical  kinetics  algorithms  to  accurately  capture 
flow  combustion  interactions  when  it  is  dominated  by  radical  concentration  and  diffusion; 

(3)  To  develop  subscale  fundamentally  based  models  which  can  be  used  to  obtain  large 
eddy  simulations  using  vortex  methods.  Subscale  models,  in  this  case,  will  be  obtained  through 
the  application  of  the  renormalization  group  theory  to  the  equations  governing  vorticity 
stretching  and  tilting  on  scales  lower  than  numerically  resolvable  scales.  The  application  of 
RNG  in  physical  vorticity  space  is  motivated  by  our  solutions  which  show  that  the  properties  of 
vortex  lines  at  the  small  scales  follow  closely  the  RNG  predictions. 

(4)  To  modify  the  transport  element  method;  a  Lagrangian  scheme  which  we  developed 
to  simulate  reacting  species  transport,  to  act  as  a  “coupling”  algorithm  between  the  solutions 
obtained  for  the  outer  flow  and  the  inner  flow. 

The  developed  methodology  will  be  applied  to  study  mixing  and  combustion  in  reacting 
shear  flow  at  high  Reynolds  numbers.  Effect  of  strong  density  variations,  high  heat  release  rates 
at  elevated  Damkohler  numbers,  and  high  Mach  numbers  will  also  be  investigated. 
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II.  PERSONNEL 

During  the  period  of  1993-1994,  funding  was  used  to  support  the  work  of  the  following 
students: 

(1)  Marios  Soteriou,  partially  supported  by  funds  from  this  contract  as  a  post  doctor. 

(2)  Van  Luu,  a  Ph.D.  student  working  on  distributed  reaction  zone  combustion  model. 

(3)  Constantin  Petrov  ,  a  Ph.D.  student  working  on  the  thin  flame  detailed  chemistry  and 
diffusion  combustion  model. 

in.  WORK  STATUS 

Overview: 

The  objectives  of  this  program  are  to  develop  a  reaction  zone  model,  or  several  reaction 
zone  models  applicable  under  different  turbulent  combustion  interaction  regimes,  and  implement 
one  or  several  of  these  models  within  a  framework  which  is  aimed  at  the  accurate  and  efficient 
simulation  of  turbulent  combustion.  Since  most  often,  the  results  of  these  simulations  are 
intended  to  predict  the  volumetric  burning  rate,  combustion  stability,  and  emission  of  minor 
species  (NOx,  CO  and  HC),  it  is  important  that  the  simulation  includes  a  representation  of  the 
detailed  fluid  dynamics  and  mixing,  the  dynamic  impact  of  heat  release  on  the  flow  dynamics, 
and  most  significantly  the  detailed  chemical  kinetics  and  transport  models  necessary  to  evaluate 
the  rate  of  reaction  and  the  concentration  of  minor  species.  The  plan  of  our  research  has  thus 
been  organized  around  addressing  these  issues  in  detail. 

The  flow  simulation  approach  has  been  developed  and  tested  before.  We  have 
demonstrated  that  using  vortex  methods,  we  can  capture  the  large  and  small  scale  structures  of  a 
turbulent  flow  (at  least  in  an  unconfined  shear  flow  environment).  The  development  of  a 
reaction  zone  model  for  high  Reynolds  number,  high  Damkohler  number  turbulent  combustion  is 
described  in  detail  in  Section  in.I  and  Appendix  1.  The  model  is  based  on  the  notion  that  in  this 
regime,  the  entire  flame  structure  thickness  is  less  than  or  equal  to  the  small  scales  of  turbulence 
and  thus  the  impact  of  turbulence  on  the  combustion  process  can  be  expressed  in  terms  of  an 
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unsteady  strain  rate  acting  on  the  plane  of  the  flame.  This  model,  while  treating  a  single  flame,  is 
capable  of  accounting  for  the  interaction  between  neighboring  flames  as  well.  The  model  has 
been  expanded  to  allow  the  implementation  of  a  detailed  chemical  kinetics  mechanism  and  the 
differential  diffusion  of  heat  and  different  species. 

In  the  coming  year,  we  will  focus  on  the  implementation  of  this  model  into  the  vortex 
simulation  of  the  reacting  shear  layer.  In  preparation  for  this  step,  we  performed  a  comparison 
between  a  vortex  simulation  of  the  same  flow  using  direct  simulation  of  the  combustion  process, 
i.e.,  without  imposing  any  assumptions  on  the  flame  structure  or  its  interactions  with  the  flow, 
and  a  simulation  in  which  the  chemistry  was  assuming  to  occur  infinitely  fast.  The  comparison 
showed  that  while  the  thin  flame  model  simplifies  the  computations  significantly,  it  is  still 
capable  of  capturing  the  essential  dynamics  of  heat  release  and  how  it  affects  the  flow  and 
mixing.  The  comparison  between  the  results  describing  shear  layer  dynamics  using  a  thin  flame 
model  and  those  using  direct  simulation  of  the  flame  structure  are  described  in  Sections  III.2  and 
in.3.1  and  in  more  detail  in  Appendices  II  and  HI.  Finally,  a  detailed  analysis  of  the  vorticity 
dynamics  of  a  reacting  shear  layer  with  high  rates  of  heat  release  and  finite  rate,  although  fast 
kinetics,  is  presented  in  Sections  III.3  and  Appendix  IE. 

IE.l  Development  of  an  Unsteady  Strained.  Complex  Chemistry  and  Transport  Flame  Model 

for  Turbulent  Combustion  Simulation 

We  have  developed  an  unsteady  strained  flame  model  with  detailed  chemical  kinetics  and 
molecular  transport.  The  model  can  be  used  independently  to  study  the  response  of  a  flame  to 
external  unsteady  flows  or  as  a  combustion  submodel  in  a  turbulent  combustion  simulation.  As  a 
model  of  the  reaction  zone  structure  in  a  turbulent  combustion  simulation,  it  is  applicable  under 
conditions  of  high  Reynolds  number  and  high  Damkohler  number  where  the  combined 
convection-diffusion  and  reaction  zones  is  smaller  than  or  equal  to  the  small  scales  of  the  flows. 
This  is  the  regime  that  has  commonly  been  known  as  the  flamelet  regime.  Our  approach  differs 
from  the  classical  flamelet  regime  in  two  significant  aspects: 
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(1)  we  treat  the  problem  as  being  unsteady  since  it  has  been  shown  that  the  flame  response  to 
a  time  dependent  strain  is  most  instantaneous;  and, 

(2)  the  coupling  between  the  flow  and  the  flame  structure  is  described  in  terms  of  the 
boundary  conditions  imposed  by  the  flow  across  the  flame  thickness  as  well  as  the  strain, 
i.e.,  we  do  not  consider  the  flame  as  a  reaction  zone  in  a  simple  stagnation  point  flow. 

As  described  in  detail  in  Appendix  I,  the  formulation  of  the  model  is  based  on  the  unsteady 
strained  laminar  flame  equations,  while  maintaining  the  strain  rate  constant  throughout  the  flame 
structure.  Following  a  series  of  mathematical  transformations,  the  equations  governing  the  flame 
structure  are  reduced  to  a  set  of  reaction-diffusion  equations  which  are  much  easier  to  integrate 
than  the  original  equations  using  a  distribution  of  points  that  adapt  to  the  time  dependent  flame 
structures.  The  model  accommodates  an  arbitrary  chemical  kinetic  algorithm  using  an  interface 
similar  to  that  adopted  in  CHEMKIN.  Results  show  that  the  model  predictions  agree  very  well 
with  those  obtained  from  the  detailed  solution  of  the  strained  flame  equations,  while  taking  a 
small  fraction  of  the  computer  time. 

During  the  coming  year  we  will  work  on  coupling  this  model  with  the  flow  simulation 
code,  develop  the  proper  boundary  conditions  and  interface  tracking  algorithms  necessary  to 
ensure  the  coupling,  and  use  this  approach  to  simulate  turbulent  combustion  in  a  shear  layer. 

m.2  Effect  of  Exothermicitv  and  Forcing  on  the  Dynamics  of  Reacting  Shear  Flow:  Infinite 
Chemistry  Case 

In  this  work,  the  effects  of  combustion  heat  release  on  a  spatially-developing,  high  Reynolds 
number,  non-premixed,  reacting  shear-layer  have  been  investigated  under  both  forced  and  unforced 
inlet  conditions  using  the  results  of  two-dimensional  numerical  simulations.  In  this  work,  combustion 
is  modeled  using  an  infinite  reaction  speed  model.  The  numerical  solution  is  obtained  using  the 
Lagrangian  transport-element  method.  Results  indicate  that  in  an  unforced  flow,  heat  release 
decreases  the  size  of  the  mixing  region  and  the  amount  of  product  formed  over  most  of  the  domain 
via:  (i)  a  delay  of  the  onset  of  the  flow  instability  and  (ii)  a  reduction  in  the  overall  spinning  and  an 
apparent  suppression  of  the  interactions  of  the  eddies.  With  external  forcing,  which  promotes  early 
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destabilization  of  the  layer,  similar  trends  of  mixing  and  combustion  reduction  are  experienced,  as 
manifested  by  reduced  eddy  spinning,  which  keeps  the  eddy  major  axis  aligned  with  the  flow 
direction,  and  alteration  of  the  eddy  pairing  interaction  to  tearing  of  smaller  eddies  by  their  larger 
neighbors.  Volumetric  (area)  expansion  weakens  the  vorticity  field  and  is  responsible  for  delaying 
the  onset  of  the  instability  in  the  unforced  case,  the  reduction  of  the  spinning  of  the  eddies  in  both 
cases  and  the  alignment  of  the  major  axis  of  the  downstream  larger  eddy  with  the  flow  direction  in  the 
forced  case.  The  suppression  of  eddy  pairing  is  attributed  to  the  presence  of  baroclinic  vorticity 
generation.  These  results  were  confirmed  by  a  simulation  in  which  the  infinite  chemistry  assumption 
was  relaxed  and  it  has  been  established,  as  a  result,  that  the  dynamic  effects  of  heat  release  can  be 
captured  using  a  thin  flame  model  (see  Section  III.3  and  Appendix  III).  The  detail  of  this  work  are 
described  in  Appendix  II. 

III.3  Vorticity  Dynamics  of  a  Reacting  Shear  Laver  with  Finite  Rate  Kinetics 

The  effects  of  combustion  on  the  vorticity  dynamics  of  a  low-Mach  number,  forced,  spatially- 
developing,  high  Reynolds  number,  reacting-shear  layer  have  been  investigated  using  the  results  of  a 
two-dimensional  simulation.  Here  the  chemical  reaction  is  modeled  by  single-step,  irreversible, 
Arrhenius  kinetics  with  high  Damkohler  number,  moderate  Karlovitz  number  and  significant  heat 
release.  The  numerical  solution  is  obtained  using  the  Lagrangian  transport-element  method.  Results 
indicate  that  a  fast  exothermic  reaction,  with  the  concomitant  density  variation,  modifies  the  shape, 
size,  speed  and  orientation  of  the  large-scale  vortical  structures  and  their  downstream  interactions, 
leading  to  an  overall  reduction  of  the  cross-stream  growth  of  the  mixing  region.  These  changes  are 
traced  to  the  two  primary  mechanisms  by  which  density  variation  due  to  heat  release  modifies 
vorticity:  volumetric  expansion  and  baroclinic  generation.  It  is  shown  that  the  weakening  of  the 
vorticity  due  to  volumetric  expansion  diminishes  the  cross-stream  mixing  zone  by  aligning  the  eddy 
major  axis  with  the  flow  direction.  Baroclinic  vorticity  generation,  on  the  other  hand,  is  responsible 
for  the  formation  of  a  band  of  positive  vorticity  on  the  outer  perimeter  of  the  large  eddies,  whose 
vorticity  is  predominantly  negative,  which  inhibits  entrainment  and  alters  the  eddy  interaction 
mechanism  from  pairing  of  adjacent  eddies  to  tearing  of  smaller  eddies  by  their  larger  neighbors. 
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Both  mechanisms  contribute  to  the  acceleration  of  the  eddies  in  the  streamwise  direction.  The  detail  of 
this  study  are  described  in  Appendix  III. 
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NOMENCLATURE 


p  -  the  mixture  density, 
v  =  (u,  v)-  the  total  velocity  of  fluid, 

V  =  (  U,  V )  -  the  inviscid  irrotational  velocity  component, 

=  ( uc,  vc ):-  the  combustion  generated  velocity, 

Yk  -  the  mass  fraction  of  species  k, 

Xk  -  the  mole  fraction  of  species  k, 

[Xk]-  the  molar  concentration  of  species  k, 

D/c  -the  mixture  averaged  diffusion  coefficient  of  species  k, 
coic  -  the  production  rate  of  species  k, 

Wk  -  the  molecular  weight  of  species  k, 
cp  -  the  constant  pressure  specific  heat  of  mixture,  mass  units, 
cp  k  -  the  constant  pressure  specific  heat  of  species  k,  mass  units, 
K+l  -the  total  number  of  chemical  species, 

T  -  the  thermodynamic  temperature, 
p  -  the  pressure, 

X  -  the  thermal  conductivity  of  mixture, 

Ru  -  the  universal  gas  constant, 

Wmijc  -  the  molecular  weight  of  mixture, 

Hk  -  the  enthalpy  of  species  p,  molar  units, 

Vk-  the  diffusion  velocity, 

Vcor  -  the  correction  velocity, 

Icf'  k  -  the  thermal  diffusion  coefficient  of  species  k, 

(x,  y,  t)  -  physical  coordinates, 

(y,t)-  coordinates  of  the  first  transformed  domain, 

(£,  t)  -  computational  coordinates, 
ef  t)  -  the  strain  rate, 

Le  k  -  the  Lewis  number  of  species  k, 
a-  the  thermal  diffusivity  of  mixture, 

Ai~  the  pre-exponential  constant  of  reaction  i, 

Ea, the  activation  energy  of  reaction  i, 

Pi-  the  temperature  exponent  of  reaction  i, 

qi-  the  rate-of-progress  variable  of  reaction  i, 

kf,  i,  kr<  i  -the  forward  and  reverse  constants  of  reaction  i, 

Kc>  i  -  the  equilibrium  constant  of  reaction  i,  concentration  units. 
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t 


6-  the  non-dimensional  temperature, 

ctid  -  third-body  coefficient  of  species  k  in  reaction  i, 

Vfo  -  the  stoichiometric  coefficient  of  species  k  in  reaction  i 
Patmr  pressure  of  one  atmosphere, 


SUBSCRIPTS 

u-  unbumt  mixture, 
b-  burnt  mixture, 
scl-  scale  value, 

*-  reference  value, 
k-  species  number, 
i-  chemical  reaction  number, 
n-  time  layer  number, 
mix-  mixture. 
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1.1  Formulation 


The  flow  configuration  of  the  model  is  depicted  in  figure  1 .  The  flame  is  located  in  the 
stagnation  point  flow  produced  by  two  counter-flowing  planar  jets  of  reactants  and 
products.  The  x-  and  y-axes  are  directed  parallel  and  perpendicular  to  the  flame, 
respectively.  The  flame  is,  in  general,  unsteady,  i.e.  its  location  and  burning  velocity 
change  with  time.  The  objectives  of  the  present  model  are  to  investigate  the  transient 
dynamics  of  the  flame  in  response  to  sudden  and  periodic  changes  in  the  strain  rate,  and 
to  develop  a  combustion  submodel  for  high  Damkohler  number  turbulent  combustion 
simulation. 

Applying  the  boundary  layer  approximation,  in  which  diffusion  along  the 
boundary  layer,  taken  here  to  be  the  x-axis,  is  negligibly  small,  we  can  write  the 
equations  governing  the  dynamics  of  flame  as  follows: 


dp  _  d(pu)  _  dipv)  _n 
dt  dx  dy  ~  ’ 


(1.1) 


dn 

dt 


+  pu 


dYk 

dx 


+  pv 


dYk 

dy 


-  -^(pYkVk)  +  COk  Wk,for  k=  0,  K 


(1.2) 


dT 

dx 


dr 


dr 


(1.3) 


P-PW~T’ 

““  mix 


(1.4) 


du 


du 


du  _  dp 


p_  +  pu_  +  pv_=.a, 


(1.5) 


where  p  is  the  mixture  density;  u  and  v  are  the  velocities  of  fluid  in  the  x  and  y 
directions,  respectively;  Y k  is  the  mass  fraction  of  species  k;  Dk  is  the  mixture  averaged 
diffusion  coefficient  of  species  k;  d>k  is  the  production  rate  of  species  k;  Wk  is  the 
molecular  weight  of  species  k;  cp  and  cPi  *  are  the  constant  pressure  specific  heats  of 
mixture  and  species  k;  (K+l)  is  the  total  number  of  chemical  species;  T  and  p  are  the 
thermodynamic  temperature  and  pressure,  respectively;  A  is  the  thermal  conductivity  of 
the  mixture;  Ru  is  the  universal  gas  constant;  Wmix  is  the  molecular  weight  of  the 
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mixture;  Hk  is  the  enthalpy  of  species  k  in  molar  units;  and  the  diffusion  velocity  Vk  is 
defined  as  follows: 


V*  =  - 


Djc  dXk 

Xk  dy 


+  V cor  + 


Dkh,k  i  dT 
Xk  T  dy 


,  for  k-0,K 


(1.5) 


where  Vcor  is  the  correction  velocity  introduced  in  order  to  satisfy  the  identity 

^  YkVk  =  0  ;  kj,  k  is  the  thermal  diffusion  coefficient.  In  the  following,  we  neglect  the 
k=0 

thermal  diffusion  velocity  ( the  last  term  of  (1.5))  and  the  correction  velocity.  As  it  was 
shown  in  [2],  this  simplification  does  not  significantly  affect  the  results.  In  the  above 
equations,  we  neglected  the  diffusion  terms  in  the  x-momentum  equation,  the  heat 
diffusion  terms  in  the  heat  and  mass  transport  equations,  and  the  dissipation  term  in  the 
energy  equation. 

Along  the  stagnation  streamline  x  =  0,  and  due  to  symmetry  : 


“  =  0’3T0' 

leading  to  the  following  simplifications  of  equations  (1.1-5): 


5p+p*l +*p.=o 

dt  dx  oy 


(1.6) 


p  +  pv  )  +  <kWk,k  =  0,K  (1.7) 

dt  dy  dy  \W mix  dy  ) 


P 


(1.9) 


while  x-momentum  equation  (1.5)  is  identically  equal  to  zero.  The  boundary  conditions 
for  equations  (1.6-8)  are: 

y  =  -~  ;  Yk=Yk,b>  T=Tb,  dT/  dy  =  0  (1.10) 


y  -  +°°  :  Yk  -  Yk>  u  ,  T  —  Tu 


(1.11) 
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where  the  subscripts  u  and  b  denote  the  values  on  the  unburnt  and  burnt  side, 
respectively.  The  boundary  conditions  for  the  velocity  field  will  be  discussed  later. 
Although  along  the  stagnation  streamline  u  =  0,  the  term  of  equation  (1.6)  does  not 

vanish.  Since  the  boundary  conditions  for  Y*  and  T  are  independent  of  x,  it  can  be 
demonstrated  [1]  that  the  flame  structure  near  the  stagnation  streamline  is  similar  to  that 
described  by  equations  (1.6-9).  Thus,  the  simplified  system  of  equations  (1.6-9)  is  able  to 
predict  the  structure  of  the  flame  in  the  vicinity  of  the  stagnation  streamline.  One 
more  simplification  of  the  problem  can  be  made.  The  energy  equation  (1.8)  contains  a 
term  of  the  form 


X  PYk^kCp,  k  £• 
k=0  y 


Smooke  [2]  showed  that  this  term  can  be  safely  neglected  virtually  for  all  important 
configurations  of  strained  and  unstrained  flame.  Eliminating  this  term,  we  can  rewrite 
equation  (1.11)  as  follows: 


pcp 


,  cT  ^  dT  .  d 


K 

X,  MkHk 

k=0 


(1.12) 


In  this  study  we  solve  the  system  of  equations  (1.6-12)  with  boundary  conditions  (1.10- 
1 1).  First,  we  decompose  the  velocity  field  into  two  components:  (1)  V  =  (  U,  V  )  ,  and 
(2)  Vc  =  (  uc,  vc ): 

v=V  +  i£  (1.1.3) 

where  7  is  the  total  velocity,  is  the  combustion  generated  velocity,  and  V  is  the 
externally  imposed  velocity  which  satisfies  the  following  incompressible  continuity 
equation: 


dU  dV_ 
dx  +  dy 


=  0 


(1.14) 


The  idea  behind  this  decomposition  is  as  follows.  Any  vector  v  can  be  represented  as  a 
sum  of  two  vectors.  The  vector  V  provides  the  velocity  field  far  away  from  the  flame. 
The  vector  accounts  for  the  change  in  velocity  induced  by  the  variation  in  density 
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within  the  reaction  zone.  It  should  be  noted  that,  in  contrast  to  seemingly  similar  models, 
our  model  does  not  employ  the  decomposition  of  the  velocity  field  into  a  mean  and  a 
disturbance  components  since,  in  the  reaction  zone,  both  velocity  components  might  be 
of  the  same  order  of  magnitude.  Nowhere  in  the  model  is  the  assumption  |v| » |vcl  used. 
Substituting  equation  (1.13)  into  continuity  equation  (1.6),  we  obtain: 


dp 

i+p 


du  dur 

dx  P~dx 


(1.15) 


Using  equation  (1.14)  and  assuming  that  on  the  stagnation  streamline  ducl  dx  =  0,v/c  get: 

dP  |  +  v^P-0  (1.16) 

dt  dy  dy 


Substitution  the  velocity  decomposition  in  equation  (1.13)  into  equations  (1.7)  and  (1.12) 
yields: 


+  p  (V  +vc  )  — # -  I  +  4  W*, *=(>,*  (1.17) 

F  dt  H  dy  dy  \Wmix  dy  ) 


pc£  +pc,(V+vc)f= 


(1.18) 


The  purpose  of  the  introduction  of  the  velocity  decomposition  is  to  apply  a  number  of 
transformations  which  can  substantially  simplify  the  solution  of  equations  (1.17-18).  If 
the  externally  imposed  inviscid  irrotational  velocity  field  V  is  specified  as 

V-(e(t)x,-e(t)y)  (119> 


where  e  is  the  strain  rate,  then  the  following  transformation  [3] 

y  =  exp({  e(t')  dt' )  y  =B(t)y,  t  =  t  (1.20) 

JO 

eliminates  the  V  d/  dy  operators  from  equations  (1.16-18).  The  physical  idea  behind  this 
transformation  is  to  translate  the  system  of  equations  into  a  reference  frame  locally 
moving  with  the  velocity  induced  by  the  constant  density  stagnation  point  flow  V.  The 
part  of  the  convective  operator  related  to  this  velocity  component  is  naturally  eliminated 
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in  this  reference  frame.  Consider  the  following  operator  dl  dt  +  V  dl  dy  in  equations 
(1.16-18).  Using  the  chain  rule  of  differentiation  and  equations  (1.20),  we  get: 

d__  dy  d  dt  d  _  E(t]  d_ 
dy  dy  dy  dy  ft 


d__  dy_ 
dt 


dt  dy  dt  dt  dy  dt 


Taking  these  equations  into  account,  the  operator  in  the  ( y,  t) 


domain  can  be  written  as: 


d/dt  +  Vd/dy  =  eB(t)y-^  +  -  +  (-ey)B(t) 

dy  dt 


d_=  d_ 
dy  dt 


Here  we  used  expression  (1.19)  for  the  externally  imposed  constant  density  velocity  field 
V .  Following  the  application  of  the  transformation  in  equation  (1.20),  the  governing 
equations  take  the  form: 


dP_  ,  B  %Pvc)  _  ^  (1.21) 

dt  dy 


P 


dt 


dYk  2  d  I  pDk  d(WmixYk)  \ 
dy  dy\Wmix  dy  I 


+  6)k  Wk,  k  =0,  K 


(1.22) 


pc*+pCpVcS*--B>±[^)-UHt  d.23) 

The  combustion  generated  velocity  vc  is  eliminated  by  applying  another  (Howarth) 
transformation  to  equations  (1.21-23): 


(1.24) 


where  the  overbar  denotes  the  density  normalized  by  its  value  on  the  reactants  side,  pu . 
The  transformation  implies  that: 
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%  as 


(1.25) 


(dp_ 

U 


—  d  d 
dy'  — +  — 

<?£  <9t 


(1.26) 


To  apply  this  transformation,  we  divide  the  continuity  equation  by  pu  and  integrate: 


f  ^dy'+B  p(y)  vc(y)  -  B  p(0)  vc(0)  =  0 

l  * 


(1.27) 


At  the  stagnation  point  y  =  0  the  total  velocity  must  be  zero.  This  requirement  imposes 
the  following  boundary  condition  vc( 0)  =0  and  equation  (1.27)  is  reduced  to: 


f 


^rdy’  =-B  p(y)vc(y) 
dt 


(1.28) 


Mapping  equations  (1.22-23)  onto  the  (£,  x)  domain  and  using  the  chain  rule, 
expressions  (1.25-26)  and  equation  (1.28),  we  get: 

p  =B2p—l p #WmixYk)  )  +  cokWkfork  =  0,K  (1.29) 

dx  \Wmix  d£  ) 


P 


(1.30) 


Finally,  equations  (1.29-30)  are  non-dimensionalized  using  the  length  scale  based  on  the 
thermal  diffusion  length,  i.e.  ysd=  iautsci,  where  au  is  the  thermal  diffusivity 
coefficient  of  the  mixture  on  the  reactants  side,  tSci  —  &u  !$u  >  and  Su  is  the  unstrained 
flame  burning  velocity  .  Introducing  the  non-dimensional  temperature, 
6  =  (T  -Tu)  I  (Tb  -  Tu),  we  get: 
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+  hdCOkWk  fork-0,K  ,  (1.31) 


=B2 p  !  <?/  PP*  ^WWt)  ' 
*  ®“  d$  \ 3£  ) 


dO  D2  —  /  0 
pcp^==  B2  p  -f— = 


dr 


ccu 


d$ 


Xp- 

X) 


JscL 


(Tb-Tu)  t 


K 

X  ®kHk, 


(1.32) 


k=0 


where  ^  /ysci ,  and  t  =T/tsci,  and  the  boundary  conditions  are  (1.10-11).  These 
equations  will  be  further  modified  after  we  consider  the  transport  model  in  the  next 
section. 


1.2.Transport  Model 


In  this  section  we  discuss  the  evaluation  of  the  transport  coefficients  X  and  Dk,  k  =  1,  K. 
We  used  the  S  ANDIA  Transport  [4]  and  Premixed  Flame  [5]  codes  to  study  the  variation 
of  these  parameters  with  temperature.  The  results  of  this  investigation  for  the  Lewis 
numbers  Le>k=  Xl  (p  cp  Dk  )  of  all  species  in  the  stoichiometric,  atmospheric  pressure, 
methane/  air  flame  are  presented  in  figures  2a-c.  We  employed  the  "starting"  chemical 
kinetics  mechanism  used  by  Rogg  [7].  From  the  figure  it  is  clear  that  the  Lewis  numbers 
of  all  species  are  almost  constant  with  temperature.  These  results  are  in  agreement  with 
these  obtained  e.g.  in  [2]  where  the  following  constant  values  of  the  Lewis  numbers  were 
used: 


Species  Lewis  Number 


CH4 . 0.97 

02 . 1.11 

H20 . 0.83 

C02 . 1.39 

H . 0.18 

0 . 0.70 

OH . 0.73 

H02 . 1.10 

H2 . 0.30 

CO . 1.10 

H202 . 1.12 

HCO . 1.27 

CH20 . 1.28 

CH2 . 1.00 

CH20 . 1.30 
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1.00 


l* 


n2 

The  same  constant  values  were  used  in  our  simulations. 

In  figure  3  the  thermal  conductivity  coefficient  A  is  plotted  as  a  function  of 
temperature.  The  coefficient  is  obtained  using  three  different  methods:  (1)  The 
assumption  A  p  =  const  -Xu  Pu,  (2)  the  Smooke  correlation  A  =  258  lO^Cp  (T/  TU)T 
(CGS  units),  r  =  0.7;  (3)  the  SANDIA  Transport  code  [4],  Methods  (2)  and  (3)  produce 
virtually  identical  results  over  the  entire  range  of  temperature,  while  method  (1) 
overpredicts  A  by  23  %  in  the  flame  zone.  The  reason  for  this  overprediction  can  be  seen 
in  figure  4  where  the  product  A  p  is  plotted  as  a  function  of  temperature.  In  the  flame 
zone  (  1200  -  2200  K)  the  product  is  almost  constant  and  the  approximation  works  fairly 
well,  but  in  the  lower  temperature  range  (T  <  1000  K)  the  product  increases  from  2.5  at  T 
=  1000  K  to  3.0  at  T  =  300  K.  The  only  advantage  of  A  p  =  const  -Xu  pu  approximation 
is  that  it  helps  to  reduce  the  governing  equations  to  a  particular  simple  form.  As  we  will 
see  later,  the  application  of  this  approximation  leads  to  an  overprediction  of  the  laminar 
burning  velocity  by  10  %.  If  one  wants  to  use  this  approximation,  the  best  results  would 
be  obtained  if  A  p  is  taken  equal  not  to  Xu  pu  but  to  some  intermediary  value,  say,  2.5 
g  erg  /  sec  -cm4-K.  In  the  following  we  used  the  Smooke  correlation  (2)  to  calculate  the 
thermal  conductivity  coefficient  A. 

Using  the  assumption  of  constant  Lewis  number,  Le>/fc  =  A/  (p  cpDk)  =  const, 
and  the  Smooke  correlation  A  =  258  lO^Cp  (77  Tu)r,  we  rewrite  the  diffusion  terms  in  the 
governing  equations  (1.31-32).  For  the  mass  fraction  equation: 


pDkp  = 


A  P  _  _J _ X  cP-u  P  Am 

cp  Le,  k  Pu  Le,k  Xu  CP  Pucp-u 


_1 _ A  CP'“ 

Le,k  Xu  CP 


P  Ctu  , 


(1.33) 


X 

where  we  used  the  definition  of  the  thermal  diffusivity  coefficient,  au  =  - - —  .  For 

Pu  cp,  u 

the  energy  equation: 


A  p  =  A  —  =  Xu-^~  —  rP'~-  -  ccup  Cp  u-^~ 
Pu  XuPuCP’u  Xu 


(1.34) 


Using  the  Smooke  correlation  and  the  equation  of  state,  we  evaluate: 


A  CP-  «  _ 

\T 1 

r 

A.  = 

Cp 

[rl 

Au 

[Tu\ 

y 

Xu 

CPM 

[Tu\ 

P  ,  Wmix  la 
Pu  Wmix,  u  T 


(1.35) 
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Dividing  equations  (1.31-32)  by  p  and  p  cp,  respectively,  and  using  equations  (1.35),  we 
obtain: 


dr 


B2 
Le,  k 


d_ 

% 


4=q ^-Yt> 

fit  wmix,u 


+ 


Ud  CQk  Wk 

P 


for  k  = 0 ,  K 


(1.36) 


B2 

dx 


jljL 

Cp  d$ 


W mix 
W mix,  u 


t $c  l 

P  Cp(  Tb'Tu) 


K 

X  ^Hk 

k=0 


(1.37) 


with  boundary  conditions  (1.10-1.11). 


1.3  The  Standard- State  Thermodynamic  Properties 


The  standard- state  thermodynamic  properties  of  species  k  are  given  in  terms  of  the 
polynomial  fits  to  the  specific  heats  at  constant  pressure: 


Cpjk  _  .  t  j.  .  t  2  j.  .  t3 

Ru 


=  ao,k  +  ai,kT  +  a2,kT  2  +  a3,kTJ  +  04,kT‘ 


(1.38) 


Hk  _  1  [  c 

RUT  RUT  J  p' 


kdT  '=  (1.39) 

ai,k  +  a2,kT/2  +a3,kT2/3  +  a4,kT3l  4  +  a5tkT4/ 5  +  a6tk/ T 


Sk  _  i  ( 

Ru  Ru  ]  T 

Jo 


(1.40) 


ait  k  InT  +  a2>kT  +  a3ikT2/ 2  +  a4tk  T 31 3  +  a^ik  T4!  4  +  a-jy 


where  the  thermodynamic  properties  are  given  in  molar  units.  Coefficients  a k  are  taken 
from  the  SANDIA  thermal  data  base  and  contained  in  the  "thermdat"  file.  The  format  of  a 
typical  entry  in  the  file  is  as  follows: 

CH3 

0.0284405 1E+02  0.06137974E-01  -0.0223034E-04  0.03785161E-08  -0.02452 15E- 12 
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0.16437809E+05  0.05452697E+02  0.02430442E+02  0.11124099E-01  -0.0168022E-03 
0.16218288E-07  -0.0586495E-10  0.16423781E+05  0.06789794E+02  15.04 
1.0e-00  0.0 


* 


The  name  of  a  species  is  written  on  the  first  line  of  the  entry  ;  The  next  seven  numbers 
represent  the  coefficients  aj>  k,j  =  1, ....  7  for  the  low  temperature  range  0<  T<  1000  K. 
The  following  seven  numbers  are  the  coefficients  ay  k  for  the  high  temperature  range 
1000 <  T<5000  K;  the  next  number  is  the  molecular  weight  of  the  species  in  g /  mole, 
followed  by  the  Lewis  number  of  the  species,  and,  finally,  the  thermal  diffusion 
coefficient  kj,  k-  This  structure  of  the  thermal  data  base  is  very  close  to  that  in  [6].  The 
only  difference  is  the  addition  of  the  molecular  weight,  the  Lewis  number  and  the  thermal 
diffusion  ratio.  Once  the  thermodynamic  properties  are  calculated  using  equations  (1.38- 
40),  the  properties  of  the  mixture  for  the  given  mass  fraction  vector  {Yo,  Yj,....,Yk)t  and 
temperature  T  are  evaluated  using  the  following  expressions: 


K 


CP  ~  X  Cp,  k 

k  =0 


Wmix  = 


'  K 
l£ 


0  _  P  Wmix 
H  RUT 


Yk/Wk 

(1.41) 

1 

- 1 

Yk/Wk 

(1.42) 

(1.43) 

where  cp  is  the  specific  heat  of  mixture  in  the  mass  units. 
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1.4  Chemical  Kinetics  Mechanism 


In  the  model  the  following  mechanism  of  18  elementary  reactions  and  13  reacting  species 
was  used  to  describe  the  methane  oxidation  [7]: 


SPECIES  * 

CH4  CH3  CH20  HCO  C02  CO 
H2  H  02  O  OH  H02  H20  N2  # 


CH4+H=>CH3+H2 

CH4+0H=>CH3+H20 

CH3+0=>CH2O+H 

CH20+H=>HC0+H2 

CH20+0H=>HC0+H20 

HC0+H=>C0+H2 

HCO+M=>CO+H+M 

HC0+02=>C0+H02 

C0+0H=>C02+H 

C02+H=>C0+0H 

H+02=>0H+0 

0H+Q=>H+02 

0+H2=>0H+H 

0H+H=>0+H2 

0H+H2=>H20+H 

H+H20=>OH+H2 

2*0H=>H2O+0 

H20+0=>2*0H 

H+02+M=>H02+M 

H02+M=>H+02+M 

H+0H+M=>H20+M 

H20+M=>H+0H+M 

H02+H=>2*0H 

H02+H=>H2+02 

H02+0H=>H20+02 


2.2E4 

3.0 

8527.7 

1.6E6 

2.1 

2452.2 

7.0E13 

0.0 

0.0 

2.5E13 

0.0 

40003.3 

3.0E13 

0.0 

1199.8 

2.0E14 

0.0 

0.0 

7.14E14 

0.0 

16811.7 

3.0E12 

0.0 

0.0 

4.4E6 

1.5 

-740.9 

1.6E14 

0.0 

26338.4 

1.2E17 

-0.91 

16532.0 

1.8E13 

0.0 

0.0 

1.5E7 

2.0 

7564.5 

6.7E6 

2.0 

5573.6 

1.0E8 

1.6 

3303.0 

4.6E8 

1.6 

18601.8 

1.5E9 

1.14 

O 

O 

1.5E10 

1.14 

17282.5 

2.0E18 

-0.8 

0.0 

2.77E18  -0.8 

46701.7 

2.15E22 

-2.0 

0.0 

3.87E23 

-2.0 

119397.7 

1.5E14 

0.0 

1001.4 

2.5E13 

0.0 

690.7 

2.0E13 

0.0 

0.0 

The  mechanism  is  capable  of  predicting  with  reasonable  accuracy  the  burning  velocity 
and  the  profiles  of  temperature  and  species  mass  fraction  in  lean  to  stoichiometric  flames. 
The  file  "methmech”,  which  contains  the  chemical  kinetics  mechanism,  is  read  by  an 
interpreter  function  before  the  calculations  starts.  The  format  of  the  file  is  as  follows: 

*  the  first  line  of  the  file  contains  the  word  "SPECIES"  followed  by  the  list  of 
species  participating  in  the  scheme.  The  names  of  the  species  in  the  list  are  separated  by  a 
blank  space.  The  order  in  which  the  species  appear  in  the  list  is  the  order  in  which  the 
species  data  will  be  written  in  the  chemical  kinetics  data  base.  The  list  of  the  species  is 
closed  by  a  symbol 
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*  next  lines  of  the  file  describe  the  chemical  kinetics  equations.  A  chemical 
kinetics  equation  must  contain  the  species  names  given  in  the  list.  The  left-hand  side  and 
the  right-hand  side  of  the  chemical  kinetics  equation  are  separated  by  a  symbol  "="  in  the 
case  of  a  reversible  reaction,  or  by  a  symbol  "=>"  in  the  case  of  an  irreversible  reaction. 
By  default  the  code  assumes  that  the  stoichiometric  coefficient  of  the  species  in  a  given 
reaction  is  unity;  in  the  opposite  case,  the  stoichiometric  coefficient  is  equal  to  the 
number  in  front  of  the  species  symbol.  The  chemical  kinetics  reaction  line  is  closed  by  a 
symbol 

*  the  three  numbers,  following  the  chemical  kinetics  equation,  are:  the  pre¬ 
exponential  constant  Af,  the  temperature  exponent  /?,•;  and  the  activation  energy  Ea,k  of 

reaction  i.  The  reaction  i  line  is  closed  by  a  blank  symbol. 

*  if  a  chemical  kinetics  equation  contains  symbol  "M",  it  means  that  the  reaction 
rate  is  enhanced  by  third  bodies.  The  values  of  the  third  body  efficiencies  for  particular 
species  are  specified  separately  for  each  reaction  mechanism  considered. 

Following  [6],  we  write  a  reaction  equation  in  the  general  symbolic  form: 

=  0.44) 

k=0 

where  is  the  stoichiometric  coefficient  of  species  k  in  reaction  i.  The  production  rate 
of  species  k  in  all  /  reactions  of  the  specified  reaction  mechanism  is  given  by  : 

C0k  =  ^(y'id-vki)  (L45) 

i=0 


where  <7/  is  the  rate-of-progress  variable  equal  to: 

K  K 

qi-kf.i  n™v“ -kr,ifltXklv'*‘  (146) 

k=0  k=0 


In  this  expression  [X^]  is  the  molar  concentration  of  species  k :  [X^]  —  p  Y)J  W/t  ;  the 
forward  rate  constant  k/t  /  of  reaction  i  is  calculated  as 

kf,  i  =  At  T  Pi  exp  (-  Ea,  il  Ru  T)  (1-47) 
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The  coefficients  Ai,  pi  and  Ea>  /  in  this  expression  are  provided  by  the  input  "methmech" 
file.  If  the  chemical  reaction  i  is  irreversible,  then  the  reverse  rate  constant  kr>  ,•  is  equal 
to  zero.  If  the  reaction  is  reversible,  the  equilibrium  constant  is  evaluated  first: 


v  .  _  Patm\  7- 

c'1  U«f) 


X  (v'ki-vki) 
o  exp 


(Via -Vld) 


Sk_ 

L  Ru 


Hk 


RJ 


(1.48) 


where  patm  is  the  pressure  of  1  atmosphere.  In  the  next  step,  the  reverse  rate  constant  kr>  ,• 
is  calculated: 


kr,  i  ~  kfy  i!  Kc>  i 


(1.49) 


If  a  reaction  contains  the  "third  body"  M,  the  rate-of-progress  expression  (1.46)  is 
multiplied  by 


K 

Mi  =  ^  a/d  [XjJ 


*=o 


By  default  the  code  presumes  that  a/d  of  species  k  in  reaction  i  is  equal  to  zero.  If  for 
some  species  a/d  is  not  equal  to  zero,  it  has  to  be  specified  in  the  code.  For  our  reaction 
mechanism  the  values  of  the  third  body  efficiencies  are  [7] : 

an,  i  =  1  an2,  i  =  1  ao2,  /  =  0.4  an2o,  i  =  6J  aco ,  i  =  0.75 

aco2,  i  =  1-5  och4  ,  i  =  0.54  as2,i=  0.4  (1.50) 


These  values  are  the  same  for  every  reaction  with  third  bodies. 

In  the  numerical  solution  of  equations  (1.36-37)  we  need  to  evaluate  the 
derivatives: 

d  cbjcl  d0  and  d  (by  dY j 

The  derivation  of  these  derivatives  is  a  tedious  procedure  and  we  will  give  only  the  final 
expressions: 


&-±  <*-*,> 

0IJ  i=0 


Mikfj 


dnfJ 

dY, 


Mikr,i 


dnr,  i  qi  p  OCj,  i 

dYj  Mi  Wj 


(1.51) 
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krt  i^Ir,  \M-i 


(1.52) 


dot  y  n  *  , .  kf  iTIf  [Mi  xp  * 

~\j  (VjQ-Vki)  ($i  "  j-  ^  ^mi  + 

°  z=0  m=0 


Lk 

x  «£ 


/=a 


where 


U 


tr.iK.Mi  X  fvrvii-)  Ct^J-  T;  >  ' 

1=0  R“T 


(Ok 

T 


L-  [Jf  ■  i  Vji  ^ 7 mix  V  v'.| 

^  S  "/ 


*i*.i 


vj«- 


Nt 


=  n  .{LL.  ^ mix  V  y’ 

^  r’‘b  ^  S  ' 


m 


Li 


Li  .  Li 

nril=n^/7v" 

1=0  1=0 


and  ^  V/t-  is  the  sum  of  all  stoichiometric  coefficients  on  the  right-hand- side  of  equation 
1=0 


2.1  Numerical  Solution 

The  numerical  solution  of  equations  (1.36-37)  with  boundary  condition  (1.10-11) 
is  obtained  using  a  Crank-Nicholson  finite-difference  scheme  with  an  implicit 
approximation  of  the  source  term  to  maintain  stability: 

'6?  =  {  (rHS/1*1  +  RHS?),  (2. 1 ) 

Ax  2 

where  Al  is  the  non-dimensional  time  step,  RHS  is  the  notation  for  the  right-hand-side  of 
equation  (1.37),  the  upper  index  denotes  the  time  level,  and  the  lower  index  is  the  grid 
point  number  in  (  x )  domain.  The  scheme  is  illustrated  by  the  finite-difference 
approximation  of  equation  (1.37).  The  derivations  for  equation  (1.36)  are  similar. 

At  time  level  n+1,  the  non-dimensional  source  term  of  equation  (1.37): 
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(2.2) 


_ tscl _ 

P  cp  (Tb  ~Tu) 


K 

X  ®kHk 


k=o 


is  approximated  by: 


— n+l 

A 


~n  —  n  •  n  n  •  n  n  n+j 

=  A  +  {9QjBx )(  dt  ~  Qi  +  (dQJdO)i  dd  =  A  +  {dQJdO )/  (0/ 


), 


•  _  n  ^ 

where  derivatives  ( dQJdY  )j  were  neglected  since  they  are,  as  it  was  shown  in  [1],  much 
—  n 

smaller  than  ( dQJdO )/  .  The  finite-difference  approximation  of  equation  (137)  is: 


_iJ 

Az  ~2 


^r)\+1+  Qi  +  (dQJde h  (9in+1An)  + 


H  %  h 


zL-LtcfilX+Qn 
\p  %  %  h 


(2.3) 


where  the  following  notation  was  introduced: 


cf  »  c  [TV-1  wmix._ 

c/'Cp( Tu)  WmixM 


(2.4) 


The  diffusion  term  in  equation  (2.3)  at  time  layer  n  is  approximated  as  follows: 

(%■£«*  fj= '*"*  a5) 

\  P  dt,  d£  )i  p  Zi+1-  Si-1  \  £i+l-  Si  Si  -  Si-1 ) 


where  cf  i+112  =  (  cf  /+;  +  cf  i)/  2  cf  1-7/2  =  (  cf  i  +  cf /./)/  2 

t-  n 

The  term  at  time  layer  n+l  is  approximated  in  a  similar  fashion.  The  derivative  (dQ/dd  )i 
is  equal  to 


(dQJd6)i 


K 


m,i 


k=0 


(Ok 


\n 


\  dr  p  I 
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The  finite-difference  representation  of  the  mass  fraction  equation  (1.36)  is  obtained  in  a 
similar  way.  The  finite-difference  equations  are  solved  by  a  Thomas  tridiagonal  matrix 
inversion  algorithm. 

The  grid  in  the  computational  domain  ( t,  £)  is  nonuniform.  It  is  adapted  in  such  a 
way  that  the  regions  of  high  gradients  and  large  variations  of  the  dependent  variables 
Y0,  Yj,  ...,  Yk,  0  are  continuously  resolved.  Each  time  step,  for  every  mass  fraction 
and  temperature  profile,  the  following  criteria  are  verified: 


l fi+i  -fi  \<  a  | max  (f)  -  min(f)  | 


[if!  ■{$!  P 


max  (spsD  -  min 


(2.6) 


where  a  and  (5  are  the  adaptation  tolerance  parameters  of  the  order  of  0.1;  fi  and  \df/  dqji 
are  the  values  of  a  dependent  variable  and  its  derivative  at  grid  point  i .  If  criteria  (2.6)  are 
not  satisfied  at  grid  point  i,  a  new  grid  point  is  introduced  half  way  between  the  points 
i+1  and  i.  The  values  of  mass  fractions  and  temperature  at  that  midpoint  are  calculated 
as  the  average  values: 


fi+1/2  =  (fi+l+fi)  1 2 


It  can  be  demonstrated  that  this  kind  of  interpolation  provides  the  second  order  of 
accuracy. 

At  the  onset  of  every  time  step  the  thermal  thickness  of  the  flame  is  compared 
with  the  distance  from  the  flame  to  the  upstream  and  downstream  boundaries  of  the 
computational  domain.  If  either  of  these  distances  is  smaller  than  two  flame  thicknesses, 
the  computational  domain  is  regridded  by  removing  every  second  grid  point  from  the 
mass  fraction  and  temperature  profiles  and  adding  them  to  the  boundaries  of  the  domain 
with  the  space  step  equal  to  the  mean  mesh  size  of  the  "old"  grid.  The  values  of  the  mass 
fractions  and  the  temperature  at  the  added  mesh  points  are  taken  to  be  equal  to  the 
boundary  values  of  the  "old"  mesh.  Hence,  the  zero-gradient  boundary  conditions  are 
naturally  satisfied  for  the  new  mesh.  The  total  number  of  mesh  points  is  not  changed  by 
the  regridding  procedure.  Each  time  step,  the  physical  coordinate  y-t  is  recalculated  in 
such  a  way  that  the  stagnation  point,  )>,•  =  0 ,  remains  fixed. 

The  time  step  of  the  calculations  is  the  minimum  of  the  characteristic  time  scales 
of  species  production: 
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tscl,  k  =  [Xk,  max  1 1  max  { 6)k)l  3 


(2.7) 


where  max  {tty}  is  the  maximum  of  the  production  rate  of  species  k  and  [Xki  max//  is  the 
molar  concentration  of  species  k  at  this  maximum;  and  the  computational  time  scale 
based  on  a  Courant  criterion: 


i.e. 


tscl,  Courant  —  min  (d  exp  (  -  2  Ka  x  )  ,  (2.8) 

dt  =  min  {tscl,  0’  tscl,  idscl,  2>—>tscl,  Kdscl,  Courant !  number  }  (2.9) 


where  the  number  is  specified  in  the  keyword  input  file  "inputl".  We  found  that  in  most 
of  our  runs  the  Courant  time  scale  is  the  smallest  and  actually  determines  the  time  step  of 
the  calculations.  Equation  (2.8)  shows  that  the  Courant  time  step  is  an  exponentially 
decreasing  function  of  time.  For  highly  strained  flames  (  large  Ka  ),  the  time  step 
decreases  rapidly  with  time,  exponentially  slowing  down  the  calculation.  A  way  to  solve 

this  problem  is  to  regrid  the  computational  domain  (increase  min  (^  ^]  )  whenever  the 
time  step  becomes  exceedingly  small. 


3.1  The  Structure  of  the  Code 

The  structure  of  the  Unsteady  Strained  Flame  Code  is  shown  in  Table  1.  The  code 
uses  the  freely-propagating  premixed  flame  mass  fraction  and  temperature  profiles 
generated  by  the  SANDIA  Premixed  Flame  Code  [5]  as  initial  profiles  for  the  solution  of 
the  unsteady  strained  flame  problem.  Both  the  SANDIA  code  and  our  code  employ  the 
same  chemical  kinetics  mechanism.  The  initial  conditions  are  written  in  the  file  "raw.dat" 
using  the  following  format: 

-  number  of  species  NSP  +3,  number  of  grid  points  NG,  pressure,  the  mass  flux 

-  grid  points  positions,  in  cm 

-  solution  matrix: 

Y(0, 0)  Y(l,  1) ...  Y(NSP,  0)  T(0)  flow  rate  (0) 

Y(0, 1)  Y(0, 1) ...  Y(NSP,  1)  T(l)  flow  rate  (1) 
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Y(0,  NG-1)  Y(0,  NG-1) ...  Y(NSP,  NG-1)  T(0)  flow  rate  (NG-1) 


where  Y(k,i)  denotes  the  mass  fraction  of  species  p  at  grid  point  i.  Reaction  description 
and  thermal  data  base  file  formats  were  described  in  detail  earlier  in  sections  1.4  and  1.3, 
respectively.  The  reaction  description  file  is  read  by  an  interpreter  function  "input".  This 
function  is  used  to  set  up  three  data  bases:  (1)  the  thermal  data  base  spc;  (2)  the  reaction 
data  base  reaction ;  and  (3)  the  species  data  base  elmnt_dat.  The  data  bases  are  organized 
as  global  structures.  Hence,  each  entry  in  a  data  base  can  be  "seen"  by  any  function  of  the 
code.  Three  data  input  files  "thermdat",  "methmech",  and  "raw.dat"  are  read  by  the 
functions  "therm_dat",  "input",  and  "init_prfls",  respectively.  The  functions  are  located  in 
the  library  functions  file  "flmlib2.c".  This  file  contains  the  functions  dealing  with  the 
calculation  of  the  chemical  kinetics  and  thermodynamics  related  data,  the  source  terms 
and  their  derivatives.  At  each  time  level,  for  a  given  solution  vector  (a:,  Y,  t}  ,  the 
functions  in  the  "fimlib2.c"  file  calculate: 

Hk,  $k>  Cp,  k>  Cp>  Wmix>  P>  kf,  i<  kr,  i<  Kc,  i  >  4i>  ®k 

(3.1) 

K  —  ~ 

£  co/cHk,  [XkJ,  dtOkf  do,  dcok/  dY j,  dQkl  3Yk,  901  96, 

k=0 


where  index  k  denotes  the  species  number  in  species  data  base  elmnt_dat.  The  module 
"flmlib2.c"  is  general  and  can  be  used  in  any  other  program  dealing  with  a  complex 
chemical  kinetics  mechanism.  The  pre-calculated  arrays  (3.1)  are  stored  as  the  global 
memory  class  elements  which  makes  them  "visible"  to  all  functions  of  the  code.  File 
"flmainS.c"  is  the  flame  code  itself.  It  is  supplied  by  the  pre-calculated  arrays  (3.1)  and 
"propagates"  the  mass  fractions  and  temperature  profiles  in  time.  It  contains  the 
subroutines  dealing  with: 

-  the  calculation  of  the  Thomas  tri-diagonal  inversion  coefficients, 

-  the  adaptation  of  the  grid, 

-  the  regridding  of  the  computational  domain, 

-  output, 

-  transformations  from  the  computational  to  physical  domain, 

-  the  calculation  of  the  variable  time  step, 

-  inversion. 

First,  the  flame  code  "flmainS.c"  reads  the  keyword  input  file  "inputl"  which  has  the 
following  format: 
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file  input  1 : 


1.0e+03 

2.0e-03 

1.0 

0.3 

0.4 

yes 

The  first  line  of  the  input  file  specifies  the  externally  imposed  strain  rate  in  1/sec;  second 
line  is  the  time  in  sec.  when  the  calculation  stops;  third  line  determines  the  computational 
time  step  of  the  calculation.  The  computational  time  step  is  equal  to  tsc[t  Courant  divided 
by  the  number  specified  in  the  third  line.  The  next  two  lines  give  the  adaptation  tolerance 
parameters  a  and  j3.  The  keyword  "yes"  means  that  the  calculation  will  start  not  from  the 
profiles  provided  by  the  SANDIA  Premixed  Flame  Code  [5],  but  from  the  previously 
recorded  output  of  the  Unsteady  Strained  Flame  Code.  This  is  frequently  done  to  match 
the  initial  location  of  the  flame  with  the  stagnation  point  in  order  to  compare,  for 
example,  the  responses  of  the  flame  to  different  strains.  Every  At  —  2  10~ 4  sec.  the  flame 
code  produces  the  output  files  mjrOO#  (the  major  species  and  the  temperature  profiles) 
and  mnrOO#  (the  minor  species  and  the  total  heat  release  profiles).  The  last  digit  "#"  in 
the  name  of  a  file  denotes  the  time  when  the  file  was  written  out.  Thus,  the  file  mnr003 
contains  the  minor  species  profiles  at  time  3  *  At  =  6  10'4  sec.  from  the  start  of  the 
calculations.  The  transient  integral  parameters  of  the  flame  (the  burning  velocity,  the  total 
heat  release  rate,  the  flame  location)  are  written  into  the  file  "omg.dat"  every  100  time 
steps  .  The  restart  file  "restart.dat"  is  produced  when  the  calculations  stop. 


4.1  Results  and  discussion 

The  results  obtained  using  the  described  model  are  compared  with  that  of  the 
Smooke  model  [2]  for  the  unstrained  flame;  and  the  Rogg  model  [7]  for  the  strained 
flame.  The  following  parameters  are  used  in  the  comparison: 

(a)  the  burning  velocity  Su  [cm/  sec] 

(b)  the  flame  location  yf  [cm]  defined  as  the  location  of  the  maximum  fuel 
consumption 

(c)  the  total  heat  release  rate  Q  [kJ/cmA2-sec] 
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(d)  the  shape  of  the  normalized  profiles,  where  the  normalized  distance  is  given 
by: 

y  =  pusu(y-yf)cp/ ^  (4,1) 

We  define  the  burning  velocity  Su  of  the  steady  propagating  unstrained  flame  as  follows. 
The  mass  fraction  equation  for  species  k  is  given  by: 

p  +  pv  ^-  =  for  k  =  0,  K  (4.2) 

H  dt  y  dy  dy\Wmix  dy  ) 

For  a  steady  propagating  flame  the  first  term  in  this  equation  is  equal  to  zero  and  the  mass 
flow  rate  is  conserved: 

pv  ££_=  +  i*Wk  (43) 

H  dy  dy  \Wmix  dy  f 

pv  =  puSu  =  const  (4.4) 

Using  equation  (4.4),  we  integrate  equation  (4.3)  with  respect  to  y  from  +«  to  to 
obtain: 

pv  (Yt.u-YkJ,)=  ^ 

Far  ahead  and  far  behind  the  flame,  in  the  unbumt  and  burnt  mixtures,  respectively,  the 
mixture  reaches  equilibrium,  and,  hence,  the  gradients  on  the  right-hand  side  of  equation 
(4.5)  vanish  and 

pv  (YktU-YkJ,  )=  j  (OkWkdy  (4.6) 

For  methane,  CH4,  YCh4 ,  b  =  0„  and  equation  (4.6)  yields  the  following  expression  for  the 
burning  velocity: 


23 


(4.7) 


Su  =  — TT^ - I  COCH4  WcH4dy 

pu  YcH4,  u 

This  expression  has  physical  meaning  only  when  the  flame  reaches  a  steady-state.  During 
the  transient  period,  the  presence  of  the  unsteady  term  obscures  the  interpretation  of  the 
expression.  Other  characteristics  of  the  flame  are  defined  as  follows: 


-  the  flame  location, 

-  the  total  heat  release  rate  Q  = 


K 

X  ®kHkdy 

k=0 


(4.8) 


These  two  characteristics  are  more  appropriate  for  the  description  of  the  transient 
behavior  of  a  flame. 


4.2  The  unstrained  flame 

Although  the  model  is  derived  for  the  unsteady  strained  flame,  it  can  be  used  to 
study  the  unstrained  flame  in  the  limiting  case  when  the  externally  applied  strain  rate 
tends  to  zero  or  is  much  smaller  than  the  extinction  strain.  The  strain  rate  in  our 
calculations  was  equal  to  1.0e-04  1/sec.  Two  models  were  used  to  evaluate  the  thermal 
conductivity  coefficient  A: 

Transport  model  1:  Ap  =  A uPu-  const  ,  and 

Transport  model  2:  A I  cp  -  A  (Tl  T)\  where  A  =  258  W4  g/  cm- sec  and  r  =  0.7 

In  figure  5a  the  burning  velocity  (4.7)  is  shown  as  a  function  of  time  for  the  transport 
model  1.  The  calculations  start  at  time  t  =  0  with  the  profiles  generated  by  the  SANDIA 
Premixed  Flame  Code  [5].  As  it  was  mentioned  earlier,  the  transient  values  of  the  burning 
velocity  has  little  physical  meaning  since  the  burning  velocity  (4.7)  is  calculated  using 
the  steady-state  assumption.  As  figure  5a  shows,  the  burning  velocity  starts  from  50  cm/ 
sec,  decreases,  and  in  5.e-04  sec.  reaches  a  steady-state  value  of  40.4  cm/  sec.  As  it  was 
mentioned  in  section  1.2,  Transport  Model  1  overpredicts  the  values  of  A  in  the  reaction 
zone  by  23  %.  This  is  the  reason  for  the  steady-state  burning  velocity  of  40.4  cm/sec 
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being  higher  than  the  experimental  value  of  37  cm /  sec.  Figures  5b  and  5c  show  the  total 
heat  release  rate  and  the  flame  location  as  functions  of  time.  The  pattern  of  the  total  heat 
release  rate  variation  is  similar  to  that  of  the  burning  velocity.  It  starts  from  a  value  of 
0.13  kJ /  cmA2-sec  at  time  zero  and  drops  to  a  steady-state  value  of  0.104  kJ/cmA2-sec  in 
5.e-04  sec.  The  flame  location,  yf,  after  some  short  transient  period,  starts  to  change 
linearly  with  time  as  it  is  expected  from  a  steady  propagating  premixed  flame.  The  slope 
of  the  yf(t )  line  gives  us  the  propagation  velocity  of  the  fuel  consumption  layer.  For  an 
unstrained  flame,  the  velocity  Su  is  related  to  the  velocity  of  the  fuel  consumption  layer 
by  the  mass  conservation  principle: 


r  p(  max  coch4)  dyf 

Pu  dt 


(4.9) 


Since  Transport  Model  1  does  not  predict  correctly  the  thermal  conductivity  coefficient  A 
in  the  reaction  zone,  in  the  following  simulations  we  used  only  Transport  Model  2. 

The  integral  characteristics  of  the  flame,  Su,  yf,  ht,  obtained  using  Transport 
Model  2  are  shown  in  figures  6a-c.  The  flame  starts  to  propagate  at  the  burning  velocity 
Su  =  45  cm/  sec ,  yf=  9.9  cm,ht  =  0.12  IcJ  /  cm2-sec  and  reaches  the  steady-state  values 
of  Su  =  36.6  cm/  sec  ,ht  =  0.1  kJ  I  cm2-sec  in  approximately  4.5e-04  sec.  The  flame 
location  is  a  linear  function  of  time.  The  steady-state  value  of  36.6  cm/sec  is  much  closer 
to  the  experimental  value  of  37.6  cm/sec  than  the  value  of  40.4  cm/sec  provided  by 
Transport  Model  1.  In  figures  7-15  we  compare  our  temperature  and  mole  fraction 
profiles  plotted  versus  the  normalized  coordinate  (4.1)  with  that  of  [2].  The  construction 
of  the  normalized  coordinate  (4.1)  suggests  that  a  good  agreement  in  the  shapes  of 
profiles  can  be  obtained  only  if:  (a)  the  profiles  in  the  physical  domain  are  correct,  and 
(b)  the  integral  characteristics  of  the  flame  ( Su  ,  in  particular)  are  correct. 

In  figure  7  the  temperature  profile  of  [2]  (dashed  line)  is  compared  with  the 
profile  of  the  present  model.  The  general  agreement  between  the  profiles  is  reasonably 
good,  although  the  model  slightly  overpredicts  the  temperature  in  the  post-flame  region. 
It  should  be  noted  that  the  equilibrium  values  of  the  temperature  in  the  present  model  and 
in  [2]  are  the  same  due  the  boundary  conditions  imposed.  The  observed  difference  in  the 
post-flame  temperature  can  be  explained  by  the  fact  that  the  size  of  domain  is  not  large 
enough  for  the  profile  to  reach  an  equilibrium  value.  We  can  expect  that  the  discrepancy 
in  the  temperature  profiles  will  affect  the  profiles  of  the  minor  species  in  this  region. 
Temperature  on  the  reactants  side  is  slightly  underpredicted.  In  figure  8  the  C//4  mole 
fraction  profiles  of  the  two  models  are  compared.  The  agreement  is  very  good,  although 
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our  model  produces  slightly  lower  values  for  XcHa  behind  the  reaction  zone  y  =  0  than 
these  obtained  in  Smooke's  model  [2].  This  can  be  explained  by  the  higher  values  of 
temperature  in  that  region.  In  figures  9-14  the  major  species  H2O  ,  H2  ,  O2  ,  CO,  CO2 
and  N2  mole  fraction  profiles  of  two  models  are  compared,  respectively.  The  agreement 
between  these  profiles  is  excellent  The  last  profile  presented  here  is  the  H  mole  fraction 
profile  shown  in  figure  15.  Although  the  part  of  the  profile  on  the  reactants  side  and  its 
maximum  value  of  it  are  captured  by  our  model  with  a  reasonable  accuracy,  the  post¬ 
flame  part  of  the  profile  is  lower  than  that  predicted  in  [2].  This  fact  could  be  explained 
by  the  high  sensitivity  of  H  profile  to  temperature  and  by  the  overprediction  of 
temperature  in  the  post-flame  region  by  our  model  (see  figure  7). 

We  conclude  that  the  present  model  of  the  unsteady  strained  flame  is  able  to 
predict  correctly  the  burning  velocity,  heat  release  rate,  the  profiles  of  temperature,  major 
and  minor  species  of  the  unstrained  flame  if  the  externally  imposed  strain  rate  is  small 
(e  <  100  sec1).  Slight  differences  in  the  burning  velocity  and  in  the  shapes  of  some 
profiles  (  T  and  H,  in  particular)  can  be  attributed  to  the  differences  in  the  chemical 
kinetics  mechanisms  utilized  in  [2]  and  in  our  model.  Although  in  both  model  the 
"starting"  mechanism  of  methane  oxidation  is  used,  the  reactions  of  CH4  oxidation  in  [2] 
are  written  in  a  Lindemann  fall-off  form  and  are  different  from  that  of  the  present  model. 


4.3  The  "unsteady"  strained  flame 

The  predictions  of  our  model  for  a  highly  strained  flame  of  £  =  1,000  1/sec  are 
compared  with  the  results  obtained  by  Rogg  [7].  In  both  simulations,  the  same  chemical 
kinetics  schemes  and  the  transport  models  are  used.  In  figures  16a-c  the  integral 
parameters  of  flame  are  shown  as  functions  of  time.  The  calculations  start  with  the 
SANDIA  profiles.  In  figure  16a  the  burning  velocity  Su  drops  from  44  cm/sec  to  a 
steady-state  value  of  34  cm/sec.  In  order  to  compare  the  simulations  with  different  strain 
rates,  the  flame  is  placed  initially  at  the  stagnation  point  y  =0.  Figure  16b  depicts  the 
flame  location,  yf,  as  a  function  of  time.  The  flame  starts  its  propagation  at  y  =  0.  After 
some  time  the  flame  begins  to  "feel"  the  influence  of  the  external  velocity  field  and  its 
speed  slows  down  continuously  until  it  reaches  a  steady-state  at  y  =  0.17  cm.  This  value 
is  higher  than  the  corresponding  value  of  0.1  cm  in  [7].  Figure  16c  shows  the  variation  of 
the  heat  release  rate  with  time.  It  goes  from  a  value  of  0.12  kJ/  cmA2-sec  and  reaches  a 
steady-state  value  of  0.8  kJ/cmA2-sec.  Similar  to  that  in  [1],  the  burning  velocity  and  the 
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heat  release  rate  reach  the  steady  state  much  faster  than  the  flame  location.  Next,  we 
compare  the  profiles  of  the  major  and  minor  species.  Again,  the  profiles  of  [7]  are  plotted 
broken  line  for  comparison.  Since  the  steady-state  location  of  the  flame  in  our  model  is 
different  from  that  of  [7],  the  profiles  of  [7]  were  shifted  in  such  a  way  that  the  positions 
of  the  fuel  consumption  layers  in  both  models  match.  In  figure  17  the  temperature 
profiles  are  compared.  These  exhibit  very  close  similarity,  especially  in  the  reaction  zone 
region.  In  the  post-flame  zone,  the  temperature  is  higher  in  the  Rogg’s  results  [7].  As  it 
was  mentioned  above,  the  equilibrium  values  of  the  temperatures  in  both  models  are  the 
same  by  definition  but  the  size  of  the  domain  shown  in  the  picture  is  not  large  enough  for 
our  profile  to  reach  this  value.  In  figures  18  and  19  the  CH4  and  O2  mass  fraction  profiles 
of  both  models  are  plotted.  They  are  virtually  identical.  In  figure  20  the  CO2  mass 
fraction  profiles  are  compared.  The  profiles  are  identical  in  the  reaction  zone  and  differ 
slightly  in  the  post-flame  region.  The  profile  of  CO2  predicted  in  [7]  has  slightly  higher 
values  of  CO2  mass  fraction  here.  The  same  could  be  said  about  the  CO  mass  fraction 
profiles  compared  in  figure  21.  The  profiles  match  almost  precisely  in  the  reaction  zone, 
including  the  maximum  values  of  CO  mass  fraction,  but  diverge  in  the  post-flame 
region.  The  comparison  of  the  H2O  profiles  is  presented  in  figure  22.  Again,  the  profiles 
match  almost  exactly  everywhere  in  the  field.  A  slight  discrepancy  is  observed  only  in  the 
reaction  zone.  The  previously  plotted  profiles  represent  the  temperature  and  the  major 
species.  The  model  described  in  this  paper  predicts  well  almost  all  of  them,  except  the 
CO  profile. 

A  more  severe  test  for  a  model  is  a  comparison  of  the  minor  species  profiles. 
These  species  are  present  in  tiny  concentrations.  Many  of  them  are  extremely  sensitive  to 
the  temperature  and  chemical  kinetics  mechanism.  In  figure  23  the  mass  fraction  profiles 
of  HO2  are  compared.  They  exhibit  reasonable  similarity  in  terms  of  thickness  although 
the  peak  value  is  slightly  overpredicted  by  our  model.  It  should  be  noticed  that  the  HO2 
profile  is  located  ahead  of  the  flame,  in  the  region  which  is  predicted  by  our  model 
particularly  well.  The  situation  changes  for  OH  profile  presented  in  figure  24.  While, 
again,  the  reaction  zone  side  of  the  profile  matches  well,  the  post-flame  part  of  our  profile 
is  significantly  higher  than  that  of  [7].  The  similar  trend  is  observed  for  the  H  and  O 
profiles  in  figures  25  and  26,  respectively,  while  the  H2  profile  is  predicted  reasonably 
well,  as  it  is  shown  in  figure  27.  The  similarity  among  the  trends  in  the  OH,  H  and  0 
profiles  can  be  explained  by  the  fact  that  they  are  related  to  each  other  through  the 
following  mechanism  of  H2  oxidation: 
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H2  +  0  <r>OH  +  H 


H  +  02  <r*OH  +  O 

Figures  24,  25  and  26  show  that  the  reason  why  our  model  overpredicts  the  post-flame 
zone  is  the  fact  that  the  flow  is  not  strained  as  much  as  it  should  be.  The  equilibrium 
values  of  the  mass  fractions  are  the  same  for  both  models.  A  higher  strain  would  cause  a 
shrinking  of  the  physical  domain  which,  in  turn,  would  move  OH,  O  or  H  profiles  closer 
to  the  flame,  making  them  steeper  in  the  post-flame  region.  The  reason  why  the  post¬ 
flame  region  is  understrained  could  be  the  absence  of  the  momentum  equation  in  our 
model.  As  it  was  mentioned  earlier  in  section  1.1,  the  momentum  equation  is  identically 
equal  to  zero  on  the  stagnation  streamline.  Nevertheless,  it  provides  a  boundary  condition 
which  links  the  strain  rate  on  the  reactants  and  products  sides  [7]: 

Eb  =  £u  Vpu/  Pb  (4-10> 

This  boundary  condition  follows  from  assumption  ( dpi  dx  ).»=  [dpi  dx )+«.  It  shows 
that  the  strain  rate  on  the  products  side  must  be  approximately  fpj~pb  ~2.75  times 
greater  than  that  on  the  reactants  side.  On  the  reactants  side,  and  inside  the  reaction  zone 
the  combustion  generated  velocity  component  vc  accounts  for  the  velocity  field  created 
by  the  density  variation  in  the  flame.  The  temperature,  major  and  minor  species  profiles 
are  predicted  correctly  in  that  region.  However,  far  from  the  reaction  zone  the  vc  velocity 
component  is  unable  to  provide  the  correct  asymptotic  behavior  of  the  velocity  field 
which  should  have  been  imposed  by  satisfying  the  momentum  equation.  The  physical 
domain  in  that  region  is  thus  understrained.  That  also  explains  why  the  steady-state 
location  of  the  flame  is  overpredicted  in  our  model. 

In  order  to  correct  this  flaw,  an  attempt  has  been  made  to  specify  the  right 
asymptotic  values  of  the  velocity  by  constructing  the  strain  rate  field  in  the  following 
fashion: 


e=eu,y>y* 

e  =  (£u  +  £b)  I2,y  =  y*  and  (4.11) 

£=£b,y<y* 
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Two  values  have  been  chosen:  (1)  y*  =  0  ,  and  (2)  y *  =  yf(t).  It  should  be  noted  that  the 
construction  of  the  strain  rate  field  in  this  way  does  not  violate  the  validity  of  the 
transformations  of  section  1.1.  The  computational  domain  is  split  into  two  sub-domains. 
The  values  of  all  dependent  variables,  including  the  total  velocity,  are  continuous  across 
the  boundary.  A  typical  steady-state  profile  for  the  first  choice  of  the  characteristic  y*  is 
shown  in  figure  28.  The  H  mass  fraction  profile  (  a  thick  dotted  line)  is  compared  with 
the  profile  obtained  for  the  case  of  the  uniform  strain  rate  field.  The  reaction  zone  parts  of 
the  profiles  are  virtually  identical,  while  the  post-flame  values  of  the  two- strain- zone, 
y*  =  0  profile  are  approximately  two  times  smaller  than  that  of  the  uniform  strain  profile. 
Also  shown  in  the  figure  in  a  dashed  line  the  profile  from  [7]  which  still  has  somewhat 
lower  values  in  the  post-flame  zone.  The  same  trend  was  observed  for  OH  and  O 
profiles  which  are  not  presented  here.  The  second  choice  of  y*  produces  unsatisfactory 
results.  The  steady-state  H  profile  is  presented  in  figure  29.  Now  the  profile  is 
"overstrained"  which  manifests  itself  by  a  lower  peak  value  and  a  more  "squeezed  shape. 
In  the  figure  the  profile  is  compared  with  that  for  the  uniform  strain  case  (dashed  line). 
The  same  trend  is  observed  for  other  minor  species  profiles  and  CO.  The  major  species 
(except  CO)  and  temperature  profiles  are  reasonably  well  predicted.  The  overstrained 
character  of  the  flame  is  also  demonstrated  in  figure  30  where  the  flame  location  is 
plotted  as  a  function  of  time.  The  steady  state  flame  location  of  0.02  cm  is  much  lower 
that  the  value  of  0. 1  cm  of  [7]. 


In  the  case  of  a  highly  strained  flame  we  conclude  that: 

(1)  Our  model  is  capable  of  predicting  correctly  the  profiles  of  temperature  and 
major  species  mass  fractions  in  the  reaction  zone  and  the  post-flame  region.  The  steady- 
state  flame  location  is  overpredicted. 

(2)  The  model  predicts  correctly  the  profiles  of  minor  species  in  the  reaction  zone. 
The  profiles  of  the  radicals  OH,  O,  H  are  significantly  overpredicted  in  the  post-flame 
region. 

(3)  The  predictions  of  the  model  for  the  radicals  profiles  in  the  post-flame  zone 
can  be  improved  by  the  introduction  of  a  two-zone  strain  rate  field  which  provides  the 
correct  asymptotic  velocity  values. 
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Figure  1.  Flame  propagation  in  a  stagnation-point  flow 
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Figure  2a.  The  detailed  transport-”starting”  chemical  kinetics  Lewis  numbers  of 
CO^CHiO  ,  CH4  ,  C//3  ,  HCO  as  functions  of  the  temperature  for  an  atmospheric 
pressure,  stoichiometric,  premixed  unstrained  methane-air  flame.  SANDIA  Transport 
Code. 


Figure  2c.The  detailed  transport-”starting”  chemical  kinetics  Lewis  numbers  of 
C>2,CO ,  O  ,  Hi ,  H  as  functions  of  the  temperature  for  an  atmospheric  pressure, 
stoichiometric,  premixed  unstrained  methane-air  flame.  SANDIA  Transport  Code.  Right- 
axis  values  are  the  results  of  Smooke  [2]. 
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Figure  2d.  The  detailed  transport-”starting”  chemical  kinetics  Lewis  numbers  of 
N2MO2  M2 0  ,  OH  as  functions  of  the  temperature  for  an  atmospheric  pressure, 
stoichiometric,  premixed  unstrained  methane-air  flame.  SANDIA  Transport  Code.  Right- 
axis  values  are  the  results  of  Smooke  [2]. 


Figure  3.  A  comparison  between  the  detailed  transport-"starting"  chemical  kinetics  value 
of  the  heat  conductivity  coefficient  A,  the  value  obtained  from  the  Smooke  correlation 

A  =  258  10~*cp  (T/  TU)T,  r=0.7,  and  from  the  assumption  A  p  =  const  =  Au  pu  for  an 
atmospheric  pressure,  stoichiometric,  premixed  unstrained  methane-air  flame. 
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Figure  4.  The  product  X  p  as  a  function  of  the  temperature  for  an  atmospheric  pressure, 
stoichiometric,  premixed  unstrained  methane-air  flame.  SANDIA  Transport  Code  results. 
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Figure  5a.  The  laminar  burning  velocity  of  a  freely-propagating,  premixed,  atmospheric 
pressure,  stoichiometric  flame  as  a  function  of  time.  Transport  model  A  p  =  const  =XU  pu, 
the  “starting”  chemical  kinetics  mechanism.  Strain  rate  e  =  10  s 
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Figure  5b.  The  heat  release  rate  of  a  freely-propagating,  premixed,  atmospheric  pressure, 
stoichiometric  flame  as  a  function  of  time.  Transport  model  A  p  =  const  =Ay  pu,  the 
“starting”  chemical  kinetics  mechanism.  Strain  rate  e  =  10  s'1 
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Figure  5c.  The  flame  location  of  a  freely-propagating,  premixed,  atmospheric  pressure, 
stoichiometric  flame  as  a  function  of  time.  Transport  model  Ap  =  const  =  Xu  pu,  the 
“starting”  chemical  kinetics  mechanism.  Strain  rate  £  =  10  s1 
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Figure  6a.  The  laminar  burning  velocity  of  a  freely-propagating,  premixed,  atmospheric 
pressure,  stoichiometric  flame  as  a  function  of  time.  Transport  model 

A  =  2 .55  10  cp  (T/  Tu)r,  the  “starting”  chemical  kinetics  mechanism.  Strain  rate 

£  =  10  S1 
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Figure  6b.  The  heat  release  rate  of  a  freely-propagating,  premixed,  atmospheric  pressure, 
stoichiometric  flame  as  a  function  of  time.  Transport  model  X  =  2.58  10^cp  (T/  Tu)r, 
r=0.7,  the  “starting”  chemical  kinetics  mechanism.  Strain  rate  e  =  10  sec'1. 
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Figure  6c.  The  flame  location  of  a  freely-propagating,  premixed,  atmospheric  pressure, 
stoichiometric  flame  as  a  function  of  time.  Transport  model  X  -  2.58  10^cp  (T/Tu)r, 
r=0.7,  the  “starting”  chemical  kinetics  mechanism.  Strain  rate  e  =  10  s'1 
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Figure  7.  A  comparison  between  the  computed  temperature  profiles  as  a  function  of  the 
normalized  distance  for  an  atmospheric  pressure,  stoichiometric,  methane-air  flame 

employing  a  X  =  2.58  lO^Cp  (77  Tu)r  transport  model  and  the  "starting"  chemical  kinetics 
mechanism  (solid  line);  and  the  profiles  obtained  by  Smooke  [2]  (dashed  line). 
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Figure  8.  A  comparison  between  the  computed  C//4  mole  fraction  profiles  as  a  function 
of  the  normalized  distance  for  an  atmospheric  pressure,  stoichiometric,  methane-air  flame 

employing  a  X- 2.58  lO^Cp  (77  Tu)r  transport  model  and  the  "starting"  chemical  kinetics 
mechanism  (solid  line);  and  the  profiles  obtained  by  Smooke  [2]  (dashed  line). 
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Figure  9.  A  comparison  between  the  computed  H2O  mole  fraction  profiles  as  a  function 
of  the  normalized  distance  for  an  atmospheric  pressure,  stoichiometric,  methane-air  flame 

employing  a  X  =  2.58  lO^Cp  (T/  Tu)r  transport  model  and  the  "starting"  chemical  kinetics 
mechanism  (solid  line);  and  the  profiles  obtained  by  Smooke  [2]  (dashed  line). 


Figure  10.  A  comparison  between  the  computed  H2  mole  fraction  profiles  as  a  function  of 
the  normalized  distance  for  an  atmospheric  pressure,  stoichiometric,  methane-air  flame 

employing  a  X  =  2  58  lO^Cp  (T/  Tu)r  transport  model  and  the  "starting"  chemical  kinetics 
mechanism  (solid  line);  and  the  profiles  obtained  by  Smooke  [2]  (dashed  line). 
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Figure  1 1.  A  comparison  between  the  computed  O2  mole  fraction  profiles  as  a  function  of 
the  normalized  distance  for  an  atmospheric  pressure,  stoichiometric,  methane-air  flame 

employing  a  X  =  2.58  lO^Cp  (77  Tu)r  transport  model  and  the  "starting"  chemical  kinetics 
mechanism  (solid  line);  and  the  profiles  obtained  by  Smooke  [2]  (dashed  line). 


Figure  12.  A  comparison  between  the  computed  CO  mole  fraction  profiles  as  a  function 
of  the  normalized  distance  for  an  atmospheric  pressure,  stoichiometric,  methane-air  flame 

employing  a  A  =  2.58  lO^Cp  (77  Tu)r  transport  model  and  the  "starting"  chemical  kinetics 
mechanism  (solid  line);  and  the  profiles  obtained  by  Smooke  [2]  (dashed  line). 
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Figure  13.  A  comparison  between  the  computed  CO2  mole  fraction  profiles  as  a  function 
of  the  normalized  distance  for  an  atmospheric  pressure,  stoichiometric,  methane-air  flame 

employing  a  A  =  2.58  lO^Cp  (77  TU)T  transport  model  and  the  "starting"  chemical  kinetics 
mechanism  (solid  line);  and  the  profiles  obtained  by  Smooke  [2]  (dashed  line). 
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Figure  14.  A  comparison  between  the  computed  N2  mole  fraction  profiles  as  a  function  of 
the  normalized  distance  for  an  atmospheric  pressure,  stoichiometric,  methane-air  flame 

employing  a  A  =  2 58  lO^Cp  (T/  Tu)r  transport  model  and  the  "starting"  chemical  kinetics 
mechanism  (solid  line);  and  the  profiles  obtained  by  Smooke  [2]  (dashed  line). 
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Figure  15.  A  comparison  between  the  computed  H  mole  fraction  profiles  as  a  function  of 
the  normalized  distance  for  an  atmospheric  pressure,  stoichiometric,  methane-air  flame 

employing  a  A  =  258  lO^Cp  (T/  Tu)r  transport  model  and  the  "starting"  chemical  kinetics 
mechanism  (solid  line);  and  the  profiles  obtained  by  Smooke  [2]  (dashed  line). 
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Figure  16a.  The  laminar  burning  velocity  of  a  strained,  premixed,  atmospheric  pressure, 
stoichiometric  flame  as  a  function  of  time.  Transport  model  A  =  2.58  10^cp  (T/Tu)r, 
r=0.7,  the  “starting”  chemical  kinetics  mechanism.  Strain  rate  £  =  10  3  s'1 


Time  [s] 

Figure  16b.  The  flame  location  of  a  strained,  premixed,  .atmospheric  pressure, 
stoichiometric  flame  as  a  function  of  time.  Transport  model  A  =  2.58  lO^Cp  (T/  TU)T, 
r=0.7,  the  “starting”  chemical  kinetics  mechanism.  Strain  rate  £  =  10  3  s1 
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Figure  16c.  The  heat  release  rate  of  a  strained,  premixed,  atmospheric  pressure, 
stoichiometric  flame  as  a  function  of  time.  Transport  model  A  =  2.58  10~^Cp  (T/Tu) , 
r=0.7,  the  “starting”  chemical  kinetics  mechanism.  Strain  rate  e  =  10  2  s1 


Figure  17.  A  comparison  between  the  computed  temperature  profile  for  a  highly  strained 
atmospheric  pressure,  stoichiometric,  methane-air  flame  employing  a 

A  =  2-55  lO^Cp  (T/Tu)T  transport  model  and  the  "starting"  chemistry  mechanism  (solid 

line)  and  the  profile  of  Rogg  [7]  (dashed  line).  e=103  s‘K  The  vertical  arrow  indicates 
the  location  of  maximum  fuel-consumption  rate. 
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Figure  18.  A  comparison  between  the  computed  CH 4  mass  fraction  profile  for  a  highly 
strained  atmospheric  pressure,  stoichiometric,  methane-air  flame  employing  a 

A  =  258  lO^Cp  (T/TU)T  transport  model  and  the  "starting"  chemistry  mechanism  (solid 

line)  and  the  profile  of  Rogg  [7]  (dashed  line),  £  =  10  3  s'1.  The  vertical  arrow  indicates 
the  location  of  maximum  fuel-consumption  rate. 
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Figure  19.  A  comparison  between  the  computed  O2  mass  fraction  profile  for  a  highly 
strained  atmospheric  pressure,  stoichiometric,  methane-air  flame  employing  a 

A  =  2.58  lO^Cp  (TITuy  transport  model  and  the  "starting"  chemistry  mechanism  (solid 

line)  and  the  profile  of  Rogg  [7]  (dashed  line).  £  =  10  3  s'1.  The  vertical  arrow  indicates 
the  location  of  maximum  fuel-consumption  rate. 
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Figure  20.  A  comparison  between  the  computed  CO2  mass  fraction  profile  for  a  highly 
strained  atmospheric  pressure,  stoichiometric,  methane-air  flame  employing  a 

X  =  2.58  lO^Cp  (T/  Tu)r  transport  model  and  the  "starting"  chemistry  mechanism  (solid 

line)  and  the  profile  of  Rogg  [7]  (dashed  line),  e  =  10  3  .W.  The  vertical  arrow  indicates 
the  location  of  maximum  fuel-consumption  rate. 
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Figure  21.  A  comparison  between  the  computed  CO  mass  fraction  profile  for  a  highly 
strained  atmospheric  pressure,  stoichiometric,  methane-air  flame  employing  a 

X  =  2J8  lO^Cp  (77  Tu)r  transport  model  and  the  "starting"  chemistry  mechanism  (solid 

line)  and  the  profile  of  Rogg  [7]  (dashed  line),  e  =  10  3  s'1.  The  vertical  arrow  indicates 
the  location  of  maximum  fuel-consumption  rate. 
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Figure  22.  A  comparison  between  the  computed  H20  mass  fraction  proftie  for  a  higWy 
strained  atmospheric  pressure,  stoichiometric,  methane- air  flame  employing  a 
A  =  2.55  10^cp  (T/Tu)'  transport  model  and  the  "starting"  chemistry  mechanism  (solid 
line)  and  the  profile  of  Rogg  [7]  (dashed  line),  e  =10  3  S’1.  The  vertical  arrow  indicates 
the  location  of  maximum  fuel-consumption  rate. 
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-  2  5S  in-»e„  a/Tu)'  transport  model  and  the  "starting”  chemistry  mechanism  (solid 
ine)  and  the  profile  of  Rogg  [7]  (dashed  line),  e  =  10*  s'*.  The  vertical  arrow  indicates 
he  location  of  maximum  fuel-consumption  rate. 
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Figure  24.  A  comparison  between  the  computed  OH  mass  fraction  profile  for  a  highly 
strained  atmospheric  pressure,  stoichiometric,  methane-air  flame  employing  a 

A  =  258  lO^Cp  (TITuy  transport  model  and  the  "starting"  chemistry  mechanism  (solid 

line)  and  the  profile  of  Rogg  [7]  (dashed  line),  e  =  10  ^  s'Y  The  vertical  arrow  indicates 
the  location  of  maximum  fuel-consumption  rate. 


Figure  25.  A  comparison  between  the  computed  H  mass  fraction  profile  for  a  highly 
strained  atmospheric  pressure,  stoichiometric,  methane-air  flame  employing  a 

A  =  258  lO^Cp  (T/  Tu)r  transport  model  and  the  "starting"  chemistry  mechanism  (solid 

line)  and  the  profile  of  Rogg  [7]  (dashed  line),  e  =  10  3  s'1.  The  vertical  arrow  indicates 
the  location  of  maximum  fuel-consumption  rate. 
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Figure  26.  A  comparison  between  the  computed  H2  mass  fraction  profUe  for  a  hig  y 
Strained  atmospheric  pressure,  stoichiometnc,  methane-air  flame  employ  g 
2  58  lO^Cn  (77 Tu)r  transport  model  and  the  "starting”  chemistry  mechanism  (solid 
line)  and  the  profile  of  Rogg  [7]  (dashed  line).  £  =  10  3  The  vertical  arrow  indicates 
the  location  of  maximum  fuel-consumption  rate. 


Figure  27.  A  comparison  between  the  computed  0  mass  fraction  profile  for  a  highly 
strained  atmospheric  pressure,  stoichiometric,  methane-air  flame  employing  a 

X  =  258  lO^Cp  (T/Tu)r  transport  model  and  the  "starting"  chemistry  mechanism  (solid 

line)  and  the  profile  of  Rogg  [7]  (dashed  line).  £  =  10  3  s'*.  The  vertical  arrow  indicates 
the  location  of  maximum  fuel-consumption  rate. 
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Figure  28.  A  comparison  between  the  computed  O  mass  fraction  profile  for  a  highly 
strained  atmospheric  pressure,  stoichiometric,  methane-air  flame  employing  a 

A  =  2.58  10^cp  (T/  TJr  transport  model  and  the  "starting"  chemistry  mechanism  (solid 
line)  and  the  profile  of  Rogg  [7]  (dashed  line).  £=10  3  s'*.  The  vertical  arrow  indicates 
the  location  of  maximum  fuel-consumption  rate. 
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Figure  29.  A  comparison  among  the  computed  H  mass  fraction  profile  for  a  highly 
strained  atmospheric  pressure,  stoichiometric,  methane-air  flame  employing  a 

X  =  258  lO^Cn  (T/TU)T  transport  model,  the  "starting"  chemistry  mechanism  and  two- 
zone  external  strain  field  with  y*  =  0  (dot-solid  line),  the  profile  obtained  for  the  constant 

external  strain  (solid  line),  and  the  profile  of  Rogg  [7]  (dashed  line).  £=10  3  S'1.  The 
vertical  arrow  indicates  the  location  of  maximum  fuel-consumption  rate. 
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Figure  30.  The  flame  location  of  a  highly  strained,  premixed,  atmospheric  pressure, 

stoichiometric  flame  as  a  function  of  time.  Transport  model  A  =  2.58  10~*cp  (T/  Tu)r,  the 
“starting”  chemical  kinetics  mechanism;  two-zone  strain,  y*  =  yf(t).  Strain  rate 

£u  =10  3  S'*  • 
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ABSTRACT 

The  effects  of  combustion  heat  release  on  a  spatially- 
developing,  high  Reynolds  number,  non-premixed,  reacting 
shear-layer  are  investigated  under  both  forced  and  unforced  inlet 
conditions  using  the  results  of  two-dimensional  numerical 
simulations.  Combustion  is  modeled  using  an  infinite  reaction 
speed  model.  The  numerical  solution  is  obtained  using  the 
La  gran  gi  an  transport-element  method.  Results  indicate  that  in  an 
unforced  flow,  heat  release  decreases  the  size  of  the  mixing 
region  and  the  amount  of  product  formed  over  most  of  the 
domain  via:  (i)  a  delay  of  the  onset  of  the  flow  instability  and  (ii) 
a  reduction  in  the  overall  spinning  and  an  apparent  suppression 
of  the  interactions  of  the  eddies.  With  external  forcing,  which 
promotes  early  destabilization  of  the  layer,  similar  trends  of 
mixing  and  combustion  reduction  are  experienced,  as  manifested 
by  reduced  eddy  spinning,  which  keeps  the  eddy  major  axis 
aligned  with  the  flow  direction,  and  alteration  of  the  eddy  pairing 
interaction  to  tearing  of  smaller  eddies  by  their  larger  neighbors. 
Volumetric  (area)  expansion  weakens  the  vorticity  field  and  is 
responsible  for  delaying  the  onset  of  the  instability  in  the 
unforced  case,  the  reduction  of  the  spinning  of  the  eddies  in  both 
cases  and  the  alignment  of  the  major  axis  of  the  downstream 
larger  eddy  with  the  flow  direction  in  the  forced  case.  The 
suppression  of  eddy  pairing  is  attributed  to  the  presence  of 
baroclinic  vorticity  generation. 


I.  INTRODUCTION 

Turbulent  reacting  shear  layers  represent  a  generic  and 
relatively  simple  model  of  many  combustion  systems.  They  are 
set  up  via  the  amplification  of  the  instability  of  the  shear  region 
between  two  reacting  fluids  originally  moving  at  different 
velocities,  are  essentially  two-dimensional  at  their  early  stages, 
and  are  dominated  by  large-scale  vortical  structures  [1].  These 
structures  are  the  regions  where  reactants  mix  and  bum  and  are, 
hence,  of  fundamental  importance  to  the  combustion  process. 

The  effects  of  the  combustion-heat  release  on  a  shear 
layer  have  been  of  interest  but  have  only  been  investigated 
recently.  Experimental  studies  [1-3]  indicate  that  heat  release 
reduces  the  growth  of  the  shear  layer  and  diminishes  the  rates  of 
mixing  and  burning.  Numerical  studies,  mostly  restricted  to  the 
idealized  temporally-evolving  flow  have  been  carried  out  for 
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both  non-premixed  [4,5]  and  premixed  [6]  layers  and  have,  to 
some  extend,  been  able  to  reproduce  this  behavior.  However, 
due  to  the  idealized  nature  of  the  temporal  model,  which  is 
mainly  a  consequence  of  the  introduction  of  periodic  boundary 
conditions  in  the  streamwise  flow  direction,  the  conclusions 
drawn  from  these  studies  must  be  treated  with  caution  [7,8]. 
The  non-premixed  calculations,  which  are  more  relevant  to  this 
work  assumed  a  constant  volume  combustion  process  and  were 
earned  out  under  an  infinite  activation  temperature  model  and  at 
quite  low  Damkohler  number. 

Numerical  studies  of  spatially  evolving  reacting  shear 
layers  [e.g.  7,9,10],  on  the  other  hand,  have  been  restricted  to 
cases  where  the  effects  of  combustion  on  the  flowfield  are 
negligibly  small.  An  exception  is  presented  in  Ref.[l  1]  where  a 
large  eddy  simulation  of  an  externally  forced,  compressible 
shear  layer  with  variable  transport  properties  and  Arrhenius 
kinetics  was  performed.  However,  combustion-heat  release  and 
the  activation  temperature  where  kept  low.  Nevertheless,  this 
work  was  instrumental  in  showing  that  the  effect  of  the 
temperature-dependent  transport  properties  is  small. 

In  this  paper  the  effects  of  combustion  heat  release  on  a 
spatially  evolving  reacting  shear  layer  are  numerically 
investigated  at  high  values  of  heat  release  and  for  both  forced 
and  unforced  inlet  conditions  using  a  high  resolution  numerical 
scheme.  An  advantage  of  such  an  approach  is  that  it  enables  an 
assessment  of  the  effect  of  heat  release  on  the  initial  instability  of 
the  flow  as  well  as  on  the  dynamics  downstream.  An  infinite 
reaction  speed  model  is  used  to  describe  the  reacting  field.  This 
model  was  shown,  in  earlier  work  [8,12],  to  adequately  predict 
the  effect  of  heat  release  on  the  flowfield  when  combustion  is 
fast  compared  to  the  flow.  The  unforced  layer  results  are 
presented  first.  This  is  done  in  terms  of  both  the  instantaneous 
and  average  flow  dynamics.  The  forced  layer  results  are 
subsequently  presented  in  a  similar  fashion.  Finally, 
explanations  of  the  flow  features  are  proposed  based  on  the 
vorticity  dynamics  and  the  related  convolution  of  the  flame. 


E.  FORMULATION 

The  motion  of  a  two-dimensional,  high  Reynolds 
number,  reacting  flow  is  considered.  Compressibility  effects  are 
allowed  under  the  low  Mach  number  assumption  [4].  The 
effects  of  gravity  are  neglected.  All  species  behave  as  perfect 
gases  with  equal  molar  masses  and  constant  and  equal 
diffusivities  and  specific  heats.  Combustion  occurs  via  a  single 
step,  irreversible,  mole  preserving,  infinite  rate  reaction. 

While  the  infinite  reaction  rate  model  eliminates  some  of 
the  parameters  involved  in  the  reaction  process,  it,  nevertheless, 
retains  the  enthalpy  of  reaction,  which  is  the  parameter  of 
interest  here.  Its  use  is  justified  by  our  earlier  work  in  which  a 
comparison  of  this  combustion  model  with  the  more  complicated 
and  computationally  expensive  Arrhenius  kinetics  model  was 
presented  [8,12].  Results  indicated  that  when  combustion  is 
fast  compared  to  the  flow  and,  hence,  the  flame  is  thin  and  does 
not  experience  quenching,  the  infinite  reaction  rate  model 
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the  material  line  stretch.  The  latter  is  readily  available  due  to  the 
Lagrangian  nature  of  the  scheme. 

In  the  second  fractional  integration  step,  diffusion  effects 
are  simulated  using  a  core  expansion  scheme  which  mimics  the 
diffusion  process  by  expanding  the  element  core  size. 

The  severe  distortion  of  the  flowmap  may  lead  to 
deterioration  of  the  accuracy  of  the  solution  due  to  insufficient 
core  overlap.  In  order  to  prevent  this,  but  without  sacrificing  the 
efficiency  of  the  calculation,  elements  are  continuously 
inserted/removed  in  regions  of  high  tensile/compressive  strains, 
as  assessed  by  the  distance  between  neighboring  elements.  This 
insertion/removal  is  carried  out  according  to  local  conservation 
laws. 


IV.  GEOMETRY  AND  BOUNDARY  CONDITIONS 

The  shear  layer  evolves  in  a  two-dimensional  channel  of 
height  H  and  length  Xmax,  between  two  parallel  streams  (1  for 
top  and  2  for  bottom)  of  different  velocities  which  mix 
downstream  of  a  thin  splitter  plate.  Both  streams  are  of  the  same 
density  and  temperature  and  each  carries  a  single  reactant.  The 
top  and  bottom  walls  are  modeled  as  rigid,  slip,  impermeable 
and  adiabatic  planes.  At  the  downstream  section,  a  condition  of 
vanishing  vorticity  and  SZ  gradient  is  used  as  outflow  boundary 
condition,  and  is  applied  by  removing  the  elements  which  cross 

the  x=Xmax  plane.  At  the  inlet  section,  the  velocity  and  the  X 
profiles  are  assumed  to  be  errorfunctions  of  equal  thickness 

leading  to  Gaussian  vorticity  and  X  gradient  profiles  of  a 
common  standard  deviation,  o.  In  contrast,  y  is  assumed  to  be 
uniform  and  equal  to  unity  there. 

The  eTrorfunction  profile  for  X  implies  that  the  profiles  of 
the  two  reacting  species  are  also  errorfunctions.  Profiles  of  this 
type  are  not  unlike  those  experimentally  observed  downstream 
of  the  splitter  plate.  The  same  is  true  for  the  velocity  profile. 
The  unity  value  for  the  y  variable  implies  that  temperature 
variation  is  only  a  consequence  of  the  presence  of  combustion 
products.  It  should  be  evident  that  such  an  inlet  condition 
coupled  with  the  nature  of  the  governing  transport  equation,  eq. 
(5),  and  the  rest  of  the  boundary  conditions  described  above, 
imply  that  y  is  equal  to  unity  throughout  the  domain  at  all  times. 
This  reduces  the  computational  cost  significantly. 

These  SZ  variable  profiles  enable  two  sets  of  temperature 
and  reacting  species  inlet  profiles:  (i)  A  non-reacting  set  where 
the  temperature  is  uniform  and  no  product  exists  and  (ii)  a 
reacting  set  where  the  flame  extends  all  the  way  to  the  inlet  and 
temperature  and  product  profiles  of  the  same  thickness  as  that  of 

the  X  variable  experience  a  maximum  at  the  channel  centerline. 
In  this  work,  the  second  set  is  chosen  due  to  the  fact  that  it 
bypasses  the  ignition  problem  close  to  the  inlet. 

The  channel  height,  H,  is  used  as  the  length  scale,  o  is 
scaled  in  such  a  way  that  two  wavelengths  of  the  most  unstable 
mode  of  the  uniform  density  shear  layer,  obtained  from  linear 
stability  analysis,  fit  within  the  channel  height,  i.e.  H=2x  12.8a. 
Under  this  condition,  the  channel  is  able  to  accommodate  a 
paired  eddy.  The  velocity  field  is  scaled  with  the  top  stream 
velocity,  Ui,  whereas  the  density  and  temperature  fields  by  the 
common  to  both  streams  values  of  these  properties,  p0  and  T0. 

Initialization  of  the  calculation  is  carried  out  by  assuming 
that  the  inlet  conditions  persist  throughout  the  domain.  Hence, 
transport  elements  are  distributed  over  nine  flat  material  layers 

(lines)  lying  within  the  region  of  finite  vorticity  and  X-gradient 
defined  by  the  above  profiles.  The  elements  are  of  square  area 

of  side  h=0.0195.  The  value  of  the  core  radius  is  8=0.0234, 
i.e.  h/8=0.833. 

Finally,  in  some  of  the  numerical  simulations  external 
forcing  is  implemented  at  the  inlet.  In  the  absence  of  forcing, 
the  evolution  of  the  layer  reflects  the  amplification  of 
numerically-excited  instability  waves,  resulting  in  the  "natural"- 
unforced  layer  behavior.  The  forcing  signal  consists  of  in-phase 
components  of  the  most  unstable  mode  of  the  uniform  density 


layer  and  its  subharmonic,  both  at  an  amplitude  of  0.025.  The 
interaction  of  the  two  forcing  frequencies  gives  rise  to  eddies  of 
disparate  size  depending  on  whether  the  two  modes  act 
constructively  or  destructively.  This  type  of  forcing  is  known  to 
lead  to  eddy  pairing  in  uniform  density,  non-reacting  shear 
layers  [13].  Implementation  of  the  forcing  is  achieved  by 
displacing  elements  at  the  inlet  according  to  the  forcing  signal. 


V.  RESULTS  AND  DISCUSSION 
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Numerical  simulations  were  carried  out  for  both  forced 
and  unforced  reacting  shear  layers.  The  time  step  for  all 

calculations  was  At=0. 1  while  the  length  of  the  domain  was 
Xmax=8  and  5  for  the  unforced  and  forced  cases,  respectively. 
A  shorter  domain  was  used  for  the  forced  case  since  the 
destabilization  of  the  flow  takes  place  closer  to  the  inlet  due  to 
the  external  forcing. 


For  all  cases,  the  inlet  velocity  ratio  is  r=r^=0.5  and  the 

Ui 

Reynolds  number  (and  Peclet  number  since  Pe=Re)  is 
Re  =  12800.  The  corresponding  Reynolds  number  based  on  the 


velocity  difference  across  the  layer,  AU=Ui-U2,  and  the 


vorticity-layer  original  thickness,  d  =  2a,  is  Red  =  ^  ^  -  500. 


The  mass  stoichiometry  ratio  is  <J>=1  and  the  non-dimensional 
enthalpy  of  reaction,  referred  to  in  this  text  as  "heat  release",  is 
Qo  =  6. 

To  assess  the  dynamic  effects  of  combustion-heat 
release,  the  reacting  calculations  were  repeated  but  the  density 
was  kept  constant  and  equal  to  its  inlet  value,  in  effect 
decoupling  it  from  the  temperature.  Since,  under  the 
assumptions  of  our  model,  the  variation  in  the  density  is  the  only 
means  by  which  the  reacting  field  can  influence  the  flow,  by 
keeping  the  density  invariant  the  flow  remains  insensitive  to  the 
reaction.  Furthermore,  the  fact  that  under  the  infinite  reaction 
speed  model  the  density  has  no  direct  effect  on  the  reaction  rate 
simplifies  the  assessment  of  the  dynamic  effects  of  heat  release 
by  distinguishing  them  from  chemical  effects.  When  this  is  not 
true,  a  more  sophisticated  approach  in  decoupling  the  flow  and 
the  reaction  is  necessary  [8]. 


Y,2  The  Unfcrod  Layer 


V,2,l  Instantaneous  FIo\yfigld 

Figures  1,  2  and  3  display  an  instantaneous  comparison 
between  the  uniform  and  variable  density  unforced  layers  in 

terms  of  the  vorticity,  X  and  the  temperature,  respectively.  It 
should  be  noted  that  the  apparent  disappearance  of  the  thin 
vorticity  layer  close  to  the  inlet  is  not  a  result  of  the  numerical 
simulation  but  of  the  resolution  of  the  post-processing  software. 

The  flow,  in  both  cases,  is  dominated  by  large  scale 
coherent  vortical  structures.  Their  formation  is  a  result  of  the 
destabilization  of  the  initially  flat  vorticity  layer  existing  between 
the  two  fluid  streams,  via  the  amplification  of  the  Kelvin- 
Helmholtz  instability  and  the  subsequent  non-linear  rollup  of  the 
layer  around  itself.  Density  variation  resulting  from  the 
combustion  heat  release  significantly  modifies  the  flowfield  and 
causes  a  decrease  in  the  layer  cross-stream  growth.  This 
appears  to  be  a  consequence  of  two  effects:  (i)  a  substantial 
delay  in  the  onset  of  the  flow  instability,  and  (ii)  an  alteration  of 
the  evolution  and  interactions  of  the  eddies. 

In  the  uniform  density  case,  the  eddies  contain  only 
negative  vorticity  and  appear  to  spin  more  vigorously.  Their 
interactions  are  by  pairings  i.e.  the  amalgamation  of  eddies  into  a 
larger  coherent  structures  initiated  by  the  clockwise  spiraling  of 
these  eddies  around  one  another.  In  the  variable  density  case, 
on  the  other  hand,  positive  vorticity  appears  on  the  eddy 
outskirts  and  is  progressively  entrained  into  their  cores.  Pairings 
appear  to  be  suppressed.  Rather,  the  growth  of  the  vortical 
structures  is  manifested  by  increases  in  the  size  of  single  eddies. 
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Figure  3.  Temperature  fields  of  the  cases  of  figure  1 


and  is  plotted  versus  the  streamwise  coordinate.  Equation  (16) 
indicates  that  the  product  thickness  calculated  here  is 
representative  of  the  mass  of  product  at  a  given  channel  cross- 
section.  Results  are  shown  in  figure  7. 


approach.  The  main  reason  for  this  is  the  relative  algebraic 
simplicity  of  the  resulting  averaged  equations  in  the  former  case. 

In  Favre  averaging,  a  quantity,  is  averaged  according 
to: 


t 


ft 


Figure  6.  The  time-averaged  product  density  profiles  of  the 
uniform  density,  uniform  flow-density  and  variable 
density  unforced  cases  at  x=5. 


where  the  tilde  indicates  Favre  averaging  and  the  overbar, 

Reynolds  averaging.  The  Favre  fluctuating  quantities,  \  ,  are, 
hence,  defined  as: 

5-5-i.  os) 

Favre  averaging  of  the  equations  of  motion  for  the  flow 
considered  here  results  in  the  following  correlations  of 
fluctuating  quantities : 

puV,  pu"2  and  pv"2  .  (19) 

As  pointed  out  in  Ref.[14]  the  first  term,  the  turbulent  stress 
term,  represents  the  rate  of  transfer  of  momentum  flux.  The 
second  and  third  terms  are  indicative  of  the  flow  turbulent  kinetic 
energy  per  unit  volume  in  the  streamwise  and  cross  stream 
directions,  respectively,  and  are  representative  of  the  turbulence 
intensity.  The  total  flow  turbulent  kinetic  energy  is: 


KE  = 


—  „2  ~  ,.2 

pu  +  pv 


(20) 
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Figure  7.  Product  density  thickness  evolution  with  the  stream- 
wise  coordinate  for  the  uniform  density,  uniform 
flow-density  and  variable  density,  unforced  cases. 


As  speculated  earlier,  the  product  thickness  decreases 
throughout  the  flow  domain  for  the  variable  density  case.  This 
is  due  to  the  drop  of  the  density  within  the  mixing  region  as  well 
as  the  decrease  in  the  size  of  the  mixing  region. 

V.2.3  Fluctuating  Flow 

Modification  of  the  instantaneous  flow  in  the  presence  of 
heat  release  implies  that  the  fluctuating  nature  of  the  flow  is  also 
altered.  This  is  illustrated  by  considering  the  turbulent 
correlations  resulting  from  averaging  the  equations  of  motion. 
In  a  variable  density  field  this  is  preferably  done  using  the  Favre 
(density  weighted)  approach  rather  than  the  traditional  Reynolds 


As  noted  earlier,  due  to  heat  release,  the  density 
decreases  significantly  and  this  results  in  reduced  Favre 
turbulent  quantities.  This  does  not  necessarily  imply  that  the 
fluctuations  are  diminished.  To  clarify  this  point,  in  what 
follows  the  uniform  flow-density  case,  is  provided  for 
comparison.  Since  this  latter  case  has  the  same  flowfield  as  the 
uniform-density  case  then  the  corresponding  Reynold's 
turbulent  quantities  are  identical.  Thus,  any  difference  in  the 
Favre  turbulent  quantities  between  the  two  cases  can  only  be  a 
consequence  of  the  variable  density  and  will  thus  provide  an 
indication  of  the  extent  of  its  effect  on  the  variable  density  case. 
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Figure  8.  The  Favre-averaged  turbulent  stress  profiles  of  the 
uniform  density,  uniform  flow-density  and  variable 
density  unforced  cases  at  x=5. 


X 

Uniform  density 

A 

Uniform  flow-density 

0 

Variable  density 

6 


Figure  10.  Vorticity  fields  of  the  uniform  density  (top)  and  variable  density  (bottom)  forced  cases  at  t=10. 


Figure  11.  X  fields  of  the  cases  of  figure  10.  The  white  contour  indicates  the  locus  of  the  reaction  interface  $.=0) 


Figure  12.  Temperature  fields  of  the  cases  of  figure  10. 
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down  in  the  presence  of  density  variation  for  most  of  the  domain 
despite  the  forcing  which  bypasses  the  instability  suppression 
seen  in  the  unforced  flow.  The  drop  is  much  more  significant  in 
the  region  where  interactions  of  the  vortical  structures  take  place. 
Thus,  the  modification  of  these  interactions  as  well  as  the 
properties  of  the  resulting  larger  structures  represent  the  most 
important  mechanisms  by  which  heat  release  decreases  the 
forced  shear  layer  cross-stream  growth. 


1 

0.8 

0.6 

0.4 

0.2 

0 

0.4  0.5  0.6  0.7  0.8  0.9  1  1.1  1.2 

u 


X 

Uniform  density 

0 

Variable  density 

:c  o 

I  ...  1  1  I  ....  1  ...  ■  Lul^L 


occupies  a  larger  area.  The  main  result  of  this  should  be  that  the 
significant  drop  in  the  mixing  region  thickness  seen  above  is  not 
as  damaging  to  the  rates  of  combustion. 
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Figure  15.  The  time-averaged  product  density  profiles  of  the 
uniform  density,  uniform  flow-density  and  variable 
density  forced  cases  at  x=2.5. 


Figure  13.  The  time-averaged  velocity  profiles  of  the  uniform 
and  variable  density,  forced  cases  at  x=2.5. 


Figure  16.  Product  density  thickness  evolution  with  the 
streamwise  coordinate  for  the  uniform  density, 
uniform  flow-density  and  variable  density,  forced 
cases. 


Figure  14.  Vorticity  thickness  evolution  with  the  streamwise 
coordinate  for  the  uniform  and  variable  density, 
forced  cases. 


In  figure  15,  the  product  density  profiles  of  the  uniform 
and  variable  density  cases  together  with  the  uniform  flow- 
density  case  are  presented.  Again,  the  similarity  between  these 
and  the  product  density  profiles  of  the  unforced  flow  is 
significant.  The  uniform  flow-density  case  experiences  a 
smaller  peak  than  the  variable  density  case  indicating  that  even 
though  the  mixing  region  is  smaller  for  this  case,  combustion  is 
more  efficient.  This  is  because  of  the  earlier  observation  that  the 
instantaneous  mixing  region,  while  smaller  in  width,  actually 


This  statement  is  quantified  in  figure  16  where  the 
product  thickness  is  plotted  versus  the  streamwise  coordinate. 
The  figure  clearly  indicates  that  heat  release  decreases  the  mass 
of  product  formed  throughout  the  flow  domain  but  that  most  of 
this  effect  is  associated  with  the  decrease  in  the  density  rather 
than  the  alteration  of  the  mixing.  In  agreement  with  earlier 
speculation,  the  drop  is  not  as  substantial  as  the  drop  in  the 
vorticity  thickness  would  suggest. 

V.3.3  Fluctuating  Flow 

The  fluctuating  features  of  the  flow  are  now  illustrated 
by  the  Favre  averaged  turbulent  stress,  figure  17,  and  kinetic 
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the  termination  of  reaction  due  to  unavailability  of  reactants 
which  results  from  an  overaccumulation  of  products,  can  and 
does  take  place.)  These  processes  are  certainly  valid  at  the  initial 
stages  of  the  flow  (before  eddy  formation)  and  at  the  braids.  As 
the  layer  rolls  into  eddies,  the  ability  of  the  flow  to  enhance 
burning  is  also  dependent  on  the  relative  concentrations  of 
reactants  which  are  reduced  due  to  the  accumulation  of  products. 
During  the  early  stages  of  the  rollup  the  availability  of  reactants 
is  practically  unimpaired  (through  entrainment  of  reactants  from 
the  free  streams),  and  the  accumulation  of  products  is  small. 
Thus,  in  this  region,  the  increases  in  product  formation  can  also 
be  related  to  the  straining  of  the  reaction  interface.  Further 
downstream,  however,  where  the  eddies  mature  and  start 
interacting,  accumulation  of  products  and  the  reduction  of 
reactants  concentrations  in  the  eddy  cores  are  more  prominent 
Thus,  stretching  of  the  reaction  interface  in  this  region  does  not 
necessarily  result  in  significant  increases  in  the  amount  of 
products  formed. 

Hence,  the  above  discussion  provides  the  explanation 
why  the  reduction  in  the  convolution  of  the  reaction  interface  in 
the  variable  density  case  results  in  a  decrease  in  the  product 
thickness. 


VI.  CONCLUSIONS 

Results  of  numerical  simulations  of  both  forced  and 
unforced  two  dimensional  spatially-developing,  high  Reynolds 
number,  non-premixed  reacting  shear-layers  indicate  that 
combustion  heat  release  reduces  the  efficiency  of  mixing  and, 
hence,  of  burning.  Moreover  it  lowers  the  intensity  of 
turbulence  and  accelerates  the  flow  in  the  stream  wise  direction. 

In  the  absence  of  external  forcing,  the  mixing  region  is 
reduced  for  most  of  the  domain  via  a  delay  in  the  onset  of  the 
flow  instability,  a  reduction  of  the  overall  spinning  of  the  eddies 
and  an  apparent  suppression  of  their  downstream  interactions. 

When  external  forcing  is  applied,  the  reduction  in  mixing 
is  mainly  due  to  the  reduced  spinning  of  the  eddies  which  keeps 
their  major  axes  aligned  with  the  flow  direction,  and  the 
alteration  of  the  eddy  interaction  mechanism  from  pairing  to 
tearing  of  smaller  eddies  by  their  larger  neighbors. 

The  streamwise  flow  acceleration  is  mainly  attributed  to 
the  volumetric  (area)  expansion  related  to  the  changes  in  density. 
This  expansion,  which  results  to  a  thicker  shear  region  close  to 
the  inlet,  is  responsible  for  the  delay  in  the  onset  of  the  flow 
instability  of  the  unforced  case.  Volumetric  expansion  also 
weakens  the  vorticity.  This  leads  to  the  reduction  of  the 
spinning  of  the  eddies  of  both  forced  and  unforced  flows  and  the 
alignment  of  the  major  axis  of  the  downstream  larger  eddy  with 
the  flow  direction  in  the  forced  case.  On  the  other  hand,  the 
resistance  to  pairing  seen  in  both  cases  is  attributed  to  baroclinic 
vorticity  generation.  This  mechanism  is  responsible  for  the 
formation  of  positive  vorticity  on  the  outer  perimeter  of  the  large 
eddies,  whose  vorticity  is  predominantly  negative.  The  positive 
vorticity  inhibits  the  clockwise  rotation  involved  in  the  eddy 
pairing  interaction  and  it  also  leads  to  diminished  entrainment 
from  the  free  streams. 
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ABSTRACT 

The  effects  of  combustion  exothermicity  on  the  vorticity  dynamics  of  a  low-Mach 
number,  forced,  spatially-developing,  high  Reynolds  number,  reacting- shear  layer  are 
investigated  using  the  results  of  a  two-dimensional  numerical  simulation.  The  chemical 
reaction  is  modeled  by  single-step,  irreversible,  Arrhenius  kinetics  with  high  Damkohler 
number,  moderate  Karlovitz  number  and  significant  heat  release.  The  numerical  solution  is 
obtained  using  the  Lagrangian  transport-element  method.  Results  indicate  that  a  fast 
exothermic  reaction,  with  the  concomitant  density  variation,  modifies  the  shape,  size,  speed 
and  orientation  of  the  large-scale  vortical  structures  and  their  downstream  interactions, 
leading  to  an  overall  reduction  of  the  cross-stream  growth  of  the  mixing  region.  These 
changes  are  traced  to  the  two  primary  mechanisms  by  which  density  variation  due  to  heat 
release  modifies  vorticity:  volumetric  expansion  and  baroclinic  generation.  It  is  shown  that 
the  weakening  of  the  vorticity  due  to  volumetric  expansion  diminishes  the  cross-stream 
mixing  zone  by  aligning  the  eddy  major  axis  with  the  flow  direction.  Baroclinic  vorticity 
generation,  on  the  other  hand,  is  responsible  for  the  formation  of  a  band  of  positive 
vorticity  on  the  outer  perimeter  of  the  large  eddies,  whose  vorticity  is  predominantly 
negative,  which  inhibits  entrainment  and  alters  the  eddy  interaction  mechanism  from  pairing 
of  adjacent  eddies  to  tearing  of  smaller  eddies  by  their  larger  neighbors.  Both  mechanisms 
contribute  to  the  acceleration  of  the  eddies  in  the  streamwise  direction. 
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INTRODUCTION 


Post-transitional,  exothermic,  reacting- shear  layers  are  encountered  in  many 
combustion  systems’  The  flow,  which  arises  due  to  the  intrinsic  instability  of  the  shear 
region  between  two  reacting  streams  originally  moving  at  different  velocities  is  endowed 
with  mechanisms  which  enhance  mixing  and  combustion  rates.  Experimental  studies  [1,2] 
indicate  that  the  essentially  two-dimensional  vortical  structures,  or  eddies,  which 
characterize  the  early  stages  of  the  non-reacting  flow  persist  in  its  reacting  counterpart 
despite  the  significant  effects  of  combustion  exothermicity.  Moreover,  most  of  the 
products  are  found  within  these  structures  exemplifying  their  importance  to  the  combustion 
process  and  motivating  further  studies  into  their  properties.  The  latter  are  essentially 
controlled  by  vorticity  dynamics. 

The  effects  of  combustion  exothermicity  on  a  post-transitional  shear  layer  have  been 
the  subject  of  several  studies.  Experimental  investigations  [1,3,4]  indicate  that  heat  release 
reduces  the  shear-layer  growth  and  diminishes  the  rates  of  mixing  and  combustion. 
Numerical  studies,  mostly  restricted  to  the  idealized  temporally-evolving  flow,  have  been 
carried  out  for  both  non-premixed  [5-6]  and  premixed  [7]  layers  and  have,  to  some  extend, 
been  able  to  reproduce  this  behavior.  However,  due  to  the  idealized  nature  of  the  temporal 
model,  conclusions  drawn  from  these  studies  must  be  treated  with  caution,  e.g.,  it  is  well 
known  that  asymmetric  entrainment  can  not  be  reproduced  by  this  model  [8-12]. 
Moreover,  the  non-premixed  calculations  [5,6]  were  limited  to  a  temperature-independent 
(zero  activation  energy)  reaction  rate  and  a  low  Damkohler  number. 

On  the  other  hand,  numerical  studies  of  spatially-evolving  layers  have  been 
restricted  to  cases  where  the  effects  of  combustion  on  the  flowfield  are  negligibly  small 
[e.g.,  13,14].  An  exception  is  presented  in  Ref.[15]  where  a  large-scale  simulation  of  a 
compressible  shear  layer  with  variable  transport  properties  and  Arrhenius  kinetics  was 
performed.  However,  combustion  heat  release  was  kept  low  and  the  effect  of  the  activation 
temperature  was  diminished  by  raising  the  background  temperature.  Nevertheless,  a  useful 
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conclusion  which  can  be  drawn  from  that  work  is  that  the  effects  of  the  temperature- 
dependent  transport  properties  play  a  weak  role  in  the  dynamics. 

In  this  paper  the  vorticity  dynamics  of  a  forced,  spatially-developing,  exothermic 
shear  layer  are  investigated  while  allowing  for  temperature-dependent  Arrhenius  kinetics 
and  significant  exothermicity.  Results  are  obtained  using  the  two-dimensional  transport- 
element  simulations.  The  numerical  scheme  makes  it  possible  to  use  higher  activation 
energy  and  Damkohler  number  than  previously  attempted.  The  effects  of  combustion 
exothermicity  on  the  vorticity  dynamics  are  analyzed  in  terms  of  the  contributions  of  the 
two  major  mechanisms  of  vorticity  modification  i.e.  volumetric  expansion  and  baroclinic 
vorticity  generation. 
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THE  GOVERNING  EQUATIONS 


A  two-dimensional  reacting  shear  flow  in  which  gravitational  effects  are  negligible 

and  compressibility  effects  are  limited  to  the  low  Mach  number  regime  [5]  is  considered. 

We  assume  that  all  species  behave  as  perfect  gases  with  equal  molecular  masses  and  equal 

and  constant  mass  diffusivities  and  specific  heats.  The  momentum  and  thermal  diffusivities 

are  also  constant  and  the  Lewis  number  is  equal  to  unity.  Combustion  occurs  according  to 

*  * 

a  single-step,  irreversible,  mole-preserving  Arrhenius  reaction  i.e.  <|>  TI1+TI2 

* 

where  TU  denotes  species,  i  =  1, 2  for  the  reacting  species  and  i  =  p  for  the  product,  and  <J> 
is  the  molar  stoichiometric  ratio. 

We  use  the  Helmholtz  decomposition  to  recast  the  equations  into  their  vorticity- 
streamfunction- velocity  potential  form,  and  Shvab-Zeldovich  (SZ)  variables  to  simplify  the 
scalar-transport  equations.  The  equations  of  the  SZ  variables  are  then  transformed  into  SZ- 
gradient  equations.  Since  combustion  occurs  at  a  finite  rate,  and  thus  the  SZ  variable 
solutions  alone  cannot  yield  the  complete  reacting  field,  one  reacting  variable  is  solved  for 
in  primitive  form.  The  resulting  equations,  in  non-dimensional  form,  are: 


V26e  =  -1— ,  =  -  co,  and  V2<j>p  =  0 

p  dt 

(la,b,c) 

=  -  (V  u)cok  +  VpxV~~  +  ^-V2cok 

dt  p2 

(2) 

4U-gVu-gx«ok)  +  jLv2g 

(3) 

^  =  jLvV(*^ 

(4) 

The  velocity  vector,  u=(u,v),  is  reconstructed  as  u  =  V<j>e  +  V<pp  Vx(\]/k).  To 
dimensionalize,  we  use  the  initial  thickness  of  the  vorticity  layer,  d,  the  velocity  difference 
between  the  two  streams,  AU,  and  the  initial  temperature  and  density,  T0  and  Po» 
respectively,  x  =  (x,y)  is  a  right-handed  Cartesian  coordinate  system,  t  is  time, 
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V  _  (JL  JL|  and  -d-  =  —  +  u- V.  fa  is  the  expansion  velocity  potential,  <J>p  is  the  velocity 
[dx’dyl  dt  at 

potential  related  to  the  boundary  conditions,  \|f  is  the  streamfunction,  <o  =  Vxu  is  the 
vorticity,  where  k  is  the  unit  vector  normal  to  the  plane  of  motion,  p  is  the  pressure,  p  is 
the  density  which  is  related  to  the  temperature,  T,  via  p  T  =  1.  The  Reynolds  and  Peclet 

numbers  are  Re  =  AIM  and  Pe  =  AIM  respectively,  where  v  is  the  kinematic  viscosity 
v  oc 

and  a  is  the  thermal  diffusivity.  g  =  Vp  is  the  gradient  of  the  SZ-variable,  p  =  X  or  7 

where  X  =  Yi-<j>Y2  and  y  =  T-  — Yp.  Yi  is  the  mass-fraction  of  species  i  with 

1-Hj> 

Yj  +  Y2  +  Yp  =  1,  <J)  is  the  mass  stoichiometric  ratio,  <J>  =  <j)  Mi  being  the  molar  mass 


of  species  i,  Q0  =  is  the  normalized  enthalpy  of  reaction  with  Cp  being  the  mixture 

CpT0  ^ 

specific  heat,  w  =  Af  p2  Yi  Y2  exp(-  -?jr)  is  the  reaction  rate  where  Af  =  is  the 

frequency  factor  and  T,  =  is  the  activation  temperature,  where  E,  and  R  are  the 

RTn 


activation  energy  and  the  universal  gas  constant,  respectively. 
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NUMERICAL  SOLUTION 

We  use  the  transport-element  method  to  perform  the  numerical  simulation  (for 
detail,  see  Ref.[12}.)  In  this  Lagrangian,  grid-free  scheme,  the  vorticity,  volumetric 
expansion,  SZ-gradients  and  product  mass  fraction  are  discretized  among  a  large  number  of 
elements  of  finite  area  and  overlapping  Gaussian  cores.  The  vortical  and  expansion 
contributions  to  the  velocity  are  computed  using  convolutions  over  the  fields  of  vortex  and 
expansion  elements,  respectively.  The  SZ-vanables  are  obtained  from  the  corresponding 
discrete  gradient  fields  via  a  similar  convolution.  The  contribution  of  the  boundary 
conditions  is  obtained  using  conformal  mapping  and  the  method  of  images. 

The  evolution  of  the  flow  and  scalar  fields  is  obtained  by  integrating  the  governing 
equations,  written  for  each  element,  in  two  fractional  steps.  In  the  first  step,  which 
includes  all  processes  except  diffusion,  the  elements  are  advected  with  the  local  velocity 
vector  and  the  relevant  source  terms  are  integrated,  e.g.,  the  circulation  of  an  element  is 
modified  by  the  baroclinic  term  according  to  the  cross  product  of  the  material  acceleration 
and  the  density  gradient  The  SZ-gradient  equation  is  substantially  simplified  using 
kinematical  relations  between  the  gradient  and  the  stretch  of  the  associated  material  line. 
The  integration  of  the  product-mass  fraction  equation,  carefully  performed  to  eliminate 
conventional  sources  of  numerical  errors,  modifies  the  local  concentrations  according  to  the 
reaction  source  term.  In  the  second  fractional  step,  diffusion  effects  are  incorporated  via  a 
novel  implementation  of  the  core  expansion  scheme  in  which  after  each  core  expansion  the 
strength  of  the  elements  is  recomputed  using  the  original  core  size  [12]. 

The  severe  distortion  of  the  flowmap,  which  could  lead  to  insufficient  core  overlap 
and  deterioration  of  the  numerical  accuracy  is  accommodated  by  continuously 
inserting/removing  elements  in  regions  of  high  tensile/compressive  strains,  as  assessed  by 
the  distance  between  neighboring  elements.  This  insertion/removal  of  elements  is  earned 
out  while  satisfying  local  conservation  laws. 
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GEOMETRY  AND  BOUNDARY  CONDITIONS 

The  shear  layer  evolves  downstream  of  a  thin  plate  in  a  two-dimensional  channel  of 
height  H,  and  length  Xm„  =  5H,  between  two  parallel  streams  of  unequal  velocities  but 
same  temperature  and  density  (To.  Po)>  each  stream  carrying  a  single  reactant  The  top  and 
bottom  walls  are  modeled  as  rigid,  slip,  impermeable  and  adiabatic  planes.  At  the  exit,  a 
condition  of  vanishing  gradient  is  imposed  by  deleting  the  elements  as  they  cross  the 
boundary.  At  the  inlet  section,  (1)  errorfunction  profiles  with  equal  thicknesses  are 
assumed  for  u  and  X,  rendering  the  0)  and  VX  profiles  Gaussians,  with  a  standard  deviation 
a,  (2)  y  is  taken  to  be  unity,  (3)  a  finite  amount  of  product  is  introduced,  in  the  form  of  a 
Gaussian  profile  with  a  standard  deviation,  o,  and  a  peak  value  of  0.4  to  enforce  ignition. 
The  scaling  of  c,  where  d=2a,  with  H  is  such  that  two  wavelengths  of  the  most  unstable 
mode  of  the  uniform-density  shear  layer  fit  within  the  channel  height. 

The  calculations  are  initialized  by  assuming  that  the  inlet  conditions  persist 

throughout  the  domain.  Square  elements  of  side  h=0.0195H  are  distributed  over  nine 

horizontal  material  layers  lying  within  the  support  of  the  discretized  quantities.  The  core 

radius  is  such  that  &  =  1.2. 

n 

External  forcing  is  implemented  at  the  inlet  by  displacing  the  elements  in  the  cross¬ 
stream  direction.  The  forcing  signal  consists  of  in-phase  components  of  the  most  unstable 
mode  of  the  uniform-density  layer  and  its  subharmonic,  both  at  an  amplitude 
Af  =  As  =  0.32d.  This  forcing,  which  promotes  rollup  and  pairing  in  the  uniform  density 
flow  [16],  results  in  the  formation  of  two  eddies  of  unequal  size  per  cycle.  A  large  eddy, 
referred  to  in  this  text  as  the  L-eddy,  forms  when  the  two  components  of  the  signal  act 
constructively  and  a  small  eddy,  the  S-eddy,  when  they  act  destructively  (see  figure  1). 
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RESULTS  AND  DISCUSSION 
Parameters 

The  following  parameters  are  used.  The  inlet  velocity  ratio  is  r  =  0.5,  Re  =  500, 
T  =10,  Af  =  100,  4>=1  and  Qq  =  6,  corresponding  to  Damkohler  number, 

Da  =  ^-=  1026,  where  tdif=^-  is  the  diffusion  time  scale  and  xche  =  -r^exp(-=f)  the 

Xche  a  f  fd 

chemical  reaction  time  scale,  and  Karlovitz  number,  Ka  =  =  0.49,  where  Xgw  =  is 

^flw  AU 

the  flow  time  scale.  The  vorticity  dynamics  of  the  reacting  layer  are  investigated  by  (1) 
comparing  the  results  of  constant-density  (CD)  and  variable-density  (VD)  simulations,  (2) 
analyzing  the  distributions  of  the  vorticity  modification  terms,  and  (3)  comparing  results  of 
simulations  in  which  one  of  the  vorticity  modification  mechanisms  is  removed.  Under  the 
assumptions  of  our  model,  the  effects  of  the  exothermic  energy  are  confined  to  density 
variation. 

The  flowfteld 

Figure  1  displays  a  comparison  between  the  CD  and  the  VD  layers,  both  depicted  in 
terms  of  the  vortex  elements  and  their  velocity  vectors  measured  with  respect  to  the  inlet 
mean  velocity.  In  agreement  with  experimental  observations,  the  flow  in  both  cases  is 
dominated  by  large-scale  vortical  structures  which  form  due  to  the  destabilization  of  the 
initial  vorticity  layer  via  the  amplification  of  the  Kelvin-Helmholtz  instability,  the  rollup  of 
the  fundamental  eddies,  and  their  interactions  due  to  the  subharmonic  instability.  The 
figure  shows  that  density  variation  modifies  the  evolution  and  interactions  of  the  vortical 
structures  leading  to  an  overall  reduction  of  the  cross-stream  growth. 

In  the  VD  case,  the  vortical  structures  accelerate  in  the  streamwise  direction.  At  the 
initial  stage,  the  fundamental  eddies  are  larger,  primarily  in  the  mean  flow  direction,  i.e., 
they  are  less  rounded.  In  the  latter  stages,  the  pairing-merging  of  fundamental  eddies  into 
larger,  coherent  and  rounded  structures  observed  in  the  CD  case  is  impaired  and  replaced 
by  tearing  of  the  S-eddy  between  its  two  neighboring  L-eddies.  This  results  in 
subharmonic  structures  which  are  less  rounded  and,  in  contrast  to  the  CD  case,  tend  to 
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keep  their  major  axis  aligned  with  the  streamwise  direction.  As  a  result,  the  layer  cross¬ 
stream  growth  is  significantly  inhibited. 

The  Temperature  Field 

Since  the  effect  of  combustion  on  the  flow  is  felt  solely  through  the  density  field, 
one  should  be  able  to  explain  the  various  phenomena  described  above  by  reviewing  the 
density,  or  temperature,  field  and  the  mechanisms  by  which  this  field  influences  the  flow. 
Figure  2  presents  instantaneous,  2D  maps  of  the  temperature,  which  in  the  case  of  the  VD 
flow  is  equivalent  to  inverse  density,  while  in  the  CD  case,  is  artificially  decoupled  from 
the  flow,  for  the  two  cases  of  figure  1.  It  is  clear  that  the  region  of  temperature  variation  is 
the  same  as  the  region  where  vortical  structures  exist.  This  is  because  these  structures 
control  the  mixing  of  the  two  reactants  and  hence,  are  where  combustion  takes  place.  The 
plot  also  provides  a  clear  indication  that  significant  phenomena  occur  at  scales  smaller  than 
those  of  the  vortical  structures.  Density  variation  is  more  closely  related  to  the  rollup  of  the 
material  layer  initially  positioned  between  the  two  reactants,  which  for  <t>  =  1,  coincides 
with  the  vorticity  layer.  The  density  experiences  a  minimum  in  the  vicinity  of  this  layer 
confirming  that  this  is  where  most  burning  occurs.  Thus,  analysis  of  the  flow  and  reaction 
fields  can  best  be  accomplished  by  investigating  the  dynamics  of  the  vorticity  layer. 

Vorticity  Dynamics,  Preliminaries 

The  combustion-related  density  field  can  influence  the  dynamics  in  several  ways. 

Mass  conservation  implies  that  density  decrease  causes  volumetric  expansion  which,  since 

the  flow  is  semi-confined,  causes  a  streamwise  acceleration.  The  vorticity  equation 

indicates  that  in  a  2D,  VD  field,  the  vorticity  of  a  fluid  element  is  modified  by  three 

VpxVp 

mechanisms.  Volumetric  expansion,  -(V.u)O),  baroclimc  generation,  — — — ,  ana 

viscous  diffusion,  ^-V2co.  In  an  exothermic  field,  V  u  is  generally  positive  and 
Kc 

volumetric  expansion  reduces  the  element  vorticity  while  preserving  its  circulation,  the 
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decrease  of  the  vorticity  is  equal  and  opposite  to  the  increase  in  area,  without  changes  in 
sign. 

Baroclinic  vorticity  generation  arises  from  the  misalignment  of  the  pressure  and 
density  gradients  and  is  capable  of  producing  vorticity  of  both  signs.  Since  the  pressure 
gradient  is  mainly  a  manifestation  of  the  fluid  acceleration,  vorticity  is  formed  as 
neighboring  elements  experience  differential  accelerations  due  to  their  different  density. 
Because  of  its  directional  nature,  baroclinic  vorticity  generation  is  strongly  linked  to  the 
geometry  of  the  flow  map. 

Viscous  diffusion  acts  to  smooth  out  vorticity  variations.  However,  when  the 
Reynolds  number  is  high,  viscous  diffusion  plays  a  minor  role. 

Vorticity  Distribution 

Figure  3  presents  the  vorticity  distribution  of  the  two  cases  shown  in  figure  1,  with 
blue  and  red  indicating  positive  and  negative  vorticity,  respectively.  The  most  striking 
difference  between  the  two  cases  is  the  presence  of  positive  vorticity,  in  the  VD  case,  on 
the  outer  rims  of  the  vortical  structures  whose  vorticity  is,  primarily,  of  a  negative  sign. 
Most  of  the  negative  vorticity  has  its  origin  in  the  incoming  flow.  Since  only  the  baroclinic 
torque  is  capable  of  producing  vorticity  of  both  signs,  the  positive  vorticity  must  be  a 
manifestation  of  this  mechanism. 

Baroclinic  Vorticity  Generation 

Figure  4  (top)  shows  the  baroclinic  vorticity  production  term  at  the  same  time  step 
as  before,  with  red  and  blue  indicating  positive  and  negative  rates.  The  complexity  of  this 
figure  can  be  unraveled  by  considering  the  action  in  the  neighborhood  on  a  single  eddy 
during  rollup.  This  is  done  in  figure  5  where  schematic  diagrams  of  the  geometric 
evolution  of  the  layer  are  drawn,  and  the  corresponding  baroclinic  term  is  estimated  by 
considering  the  cross  product  of  the  density  gradient  and  the  fluid  acceleration. 

During  rollup,  the  layer  acquires  a  "Z"  configuration.  Throughout  this  process, 
fluid  from  the  braids  is  accelerated  in  the  streamwise  direction  towards  the  evolving 
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structure.  This  acceleration  changes  its  direction  as  one  crosses  the  stagnation  point, 
defined  in  the  frame  of  reference  moving  with  the  eddy,  between  the  neighboring 
structures.  The  density  field  is  such  that  the  density  gradient  is  approximately  normal  to  the 
layer,  changing  sign  across  it,  while  the  local  stretch  ensures  that  the  magnitude  of  the 
gradient  stays  large.  Thus,  in  the  braids,  the  enhanced  density  gradient  and  the  fluid 
acceleration  are  mutually  normal  and  the  contribution  of  the  baroclinic  term  is  maximum. 
The  arrangement  of  the  vectors  involved  creates  positive  and  negative  vorticity  on  the  upper 
and  lower  sides  of  the  braid  on  the  left-hand  side  of  the  eddy,  respectively,  while  the 
opposite  occurs  on  the  right-hand  side. 

Closer  to  the  center  of  the  structure,  which  in  the  moving  reference  frame  is  also  a 
stagnation  point,  the  acceleration  on  both  sides  reverses  sign  and  so  does  the  generated 
vorticity.  The  baroclinic  term  becomes  weaker  as  the  fluid  acceleration  is  increasingly 
dominated  by  rotation  and  the  density  gradient  zone  thickens.  As  rollup  continues,  a  region 
of  negative  vorticity  is  established  inside  the  eddy  core,  surrounded  by  a  band  of  strong 
positive  vorticity  on  the  outside,  particularly  on  the  upper  left  and  lower  right  sides,  as 
shown  in  figures  4  and  5.  Adding  an  overall  negative  vorticity  to  this  picture  would  result 
in  the  vorticity  field  shown  in  figure  3. 

This  rearrangement  of  the  vorticity  has  significant  implications.  For  example, 
positive  vorticity  around  the  eddies  decreases  the  entrainment  of  irrotational  fluid  into  the 
large  structures,  since  it  weakens  the  field  created  by  the  original  negative  vorticity. 
Furthermore,  it  may  upset  the  eddy  interaction  mechanism  which  is  responsible  for  pairing 
in  the  CD  flow.  This  vorticity,  as  will  be  shown  in  detail  later,  also  induces  a  streamwise 
velocity  on  the  eddies. 

Divergence 

The  effect  of  the  expansion  term  on  the  vorticity  field  can  be  inferred  from  the 
vorticity  plots.  Figure  3b  shows  that  the  vorticity  occupies  more  area  and  that  the  negative 
vorticity  is  weaker  than  in  the  CD  case.  To  highlight  the  zone  where  volumetric  expansion 
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has  the  strongest  effect  on  the  vorticity,  we  plot  in  figure  4  (bottom)  the  distribution  of  the 
-(V.u)cd  term,  i.e.  the  rate  of  vorticity  change  due  to  divergence.  Since  divergence  is 

essentially  due  to  the  Lagrangian  change  in  density,  this  zone  can  also  be  considered  as 
where  the  flame  is  located.  In  this  regard,  the  figure  indicates  that  most  of  the  burning 
occurs  within  the  braids,  which  is  expected  since  Da  is  high  while  Ka  is  moderate. 
Evidently,  the  expansion  in  the  braids  weakens  the  vorticity  there  and  through  the  roll-up 
leads  to  the  formation  of  less  vigorous  eddies. 

Separating  the  Mechanisms 

The  above  discussion  indicates  that  the  approach  of  analyzing  the  effect  of  the 
different  mechanisms  of  vorticity  modification  using  simulations  which  involve  all  the 
terms  in  the  governing  equation  has  the  following  limitations:  (1)  it  is  not  obvious  how  the 
two  mechanisms  interact,  and  (2)  it  is  not  clear  whether  the  changes  in  the  evolution  of 
individual  eddies,  e.g.,  acceleration,  or  the  shift  of  the  eddy  interaction  mechanism  from 
pairing  to  tearing  is  a  result  of  expansion,  baroclinicity,  or  both.  To  overcome  these 
limitations,  the  calculation  for  the  VD  case  was  repeated  twice,  but  at  each  time,  one  of  the 
two  mechanisms  of  vorticity  modification  was  eliminated  in  an  attempt  to  distinguish 
between  the  phenomena  attributable  to  each  mechanism.  These  simulations  are  labeled 
baroclinicity-only  and  expansion-only  simulations. 

Figure  6  displays  the  results  of  these  calculations.  The  figure  shows  that  the 
suppression  of  pairing  and  the  onset  of  tearing  are  due  to  the  baroclinic  generation,  while 
the  alignment  of  the  major  axis  of  the  structure  with  the  streamwise  direction  is  due  to 
volumetric  expansion.  In  the  expansion-only  simulation,  figure  6  (bottom),  pairing  occurs 
in  a  way  similar  to  that  found  in  the  CD  flow.  In  this  case,  the  vorticity  layer  thickens  due 
to  volumetric  expansion,  and  rollup  is  slowed  down  due  to  the  weakened  vorticity. 

In  the  baroclinicity-only  case,  figure  6  (top),  on  the  other  hand,  the  S-eddies  are 
torn  between  their  neighboring  L-eddies,  similar  to  the  VD  case.  Positive  vorticity 
surrounding  the  vortical  structures  inhibits  the  spiraling  of  the  fundamental  eddies  around  a 
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common  center.  This  allows  enough  time  for  the  larger  L-eddies  to  gain  fluid  at  the 
expense  of  the  S-eddies.  The  above  argument  implies,  as  shown  before  [17]  that  in  the 
absence  of  external  forcing,  baroclinic  vorticity  generation  suppresses  eddy  interaction. 

Comparisons  between  the  fundamental  eddies  of  the  CD  and  baroclinicity-only  case 
indicate  that  these  eddies  are  smaller  and  more  elliptical  in  the  latter,  i.e.,  baroclinicity  is 
responsible  for  reducing  entrainment,  as  mentioned  earlier.  Figure  6  also  reveals  that  the 
acceleration  of  the  vortical  structures  in  the  streamwise  direction  is  not  only  due  to 
volumetric  expansion  but  also  due  to  baroclinicity.  This  has  been  observed  before  in  the 
study  of  the  non-reacting  VD  flow  [18].  This  acceleration  is  due  to  the  redistribution  of  the 
vorticity  inside  the  structures  and  the  loss  of  eddy  symmetry  which  arises  from  the  unequal 
tearing  of  the  S-eddy  towards  its  neighboring  L-eddies. 

Finally,  the  number  of  elements  used  to  discretize  the  evolving  material  layers 
indicates  that  baroclinic  vorticity  generation  increases  the  eddy-induced  strain  while 
volumetric  expansion  reduces  it 
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CONCLUSIONS 

Combustion  exothermicity  modifies  the  evolution  and  interactions  of  the  vortical 
structures  in  a  spatially-developing  shear  layer.  It  flattens  the  eddies  and  complicates  their 
vorticity  distribution,  increases  their  speed,  alters  their  interaction  mechanism  from  pairing 
to  tearing,  and  keeps  their  major  axis  aligned  with  the  flow  direction.  Consequently  the 
size  of  the  mixing  region  is  reduced. 

The  alteration  of  the  eddy  shape  is  a  consequence  of  the  baroclinic  torque  which 
significantly  complicates  the  eddy  vorticity  field  and, reduces  entrainment.  Vorticity 
generation  is  also  responsible  for  the  modification  of  the  eddy  infection  jnechamsm.  The 
alignment  of  the  eddy  major  axis  with  the  streamwise  direction  is  a  manifestation  of  the 
weakening  vorticity  field  due  to  volumetric  expansion.  The  acceleration  of  the  eddies  is  a 
combined  effect  of  both  mechanisms;  volumetric  expansion  because  of  confinement  and 
baroclinic  generation  because  it  induces  vorticity  asymmetry  within  each  individual  eddy. 

At  the  initial  stages,  baroclinic  generation  reduces  the  eddy  size  by  suppressing 
entrainment  while  volumetric  expansion  weakens  the  existing  vorticity.  Further 
downstream,  where  eddy  interactions  take  place,  both  mechanisms  act  to  decrease  the 
mixing  region;  baroclinic  generation  by  resisting  pairing  and  volumetric  expansion  by 
aligning  the  eddy  major  axis  with  the  flow  direction.  For  forced  shear  layers  this  latter 
effect  is  by  far  the  more  significant 
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FIGURE  CAPTIONS 


Figure  1.  The  flowfield  in  the  uniform-density,  top,  and  variable-density,  bottom,  cases  at 
t  =  64,  depicted  using  the  vortex  elements  and  their  velocity  vectors  measured 
with  respect  to  the  mean  flow.  The  top/bottom  number  of  elements  is 
5389/4034.  L  and  S  denote  examples  the  fundamental  eddies  which  form  during 
the  first  and  second  halves  of  the  subharmonic  forcing  function.  Arrows  indicate 
the  location  of  eddies  formed  at  the  same  time. 


Figure  2.  The  temperature  fields  of  the  uniform-density,  top,  and  the  variable-density  case, 
bottom,  at  t  =  64.  The  blue  to  red  color  bar  describes  the  temperature  scale 
1  <T  £  3.75. 

Figure  3.  The  vorticity  distribution  in  the  uniform-density,  top,  and  variable-density  case, 
bottom,  at  t  =  64.  The  blue  to  red  color  bar  describes  the  vorticity  scale 
-1.095  £©£0.235. 

Figure  4.  The  distributions  of  the  rates  of  vorticity  modification  by  baroclinic  generation, 
top,  and  volumetric  expansion,  bottom,  at  t  =  64.  The  blue  to  red  color  bar 
describes  the  vorticity-rate  scale  -0.02  £  cb  £  0.02 

Figure  5.  Schematic  illustration  of  the  process  of  baroclinic  vorticity  generation.  co0  is  the 
initial  vorticity  and  ebb  the  vorticity  generation  rate.  Vp  is  the  density  gradient 

and  a  the  material  acceleration. 

Figure  6.  The  baroclinicity-only,  top,  and  expansion-only,  bottom,  cases  at  t  =  64  shown 
in  terms  of  the  vortex  elements  and  their  relative  velocity.  The  top/bottom 
number  of  elements  is  9232/3818.  Arrows  indicate  the  location  of  eddies  formed 
at  the  same  time.  The  lines  on  the  left/right  sides  of  the  arrows  correspond  to  the 
locations  of  the  same  eddies  in  the  uniform/variable  density  cases,  shown  in 
figure  1. 
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